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ABSTRACT 


A  stable,  convergent  numerical  model  has  been  developed  to 
study  the  evolution  of  mesoscale  meanders  of  oceanic  fronts.  A  finite 
difference  model  was  used  with  the  frontal  zone  divided  into  two 
regions,  a  dissipative  zone  near  the  surface  expression  of  the  front, 
and  an  inviscid  region  where  the  frontal  interface  was  deeper.  \  Tests 
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using  this  model  focused  cn  some  of  the  parameters- nvhiciraf feet  the 

/ 

stability  characteristics  of  these  meanders.^The  parameters  studied 
were:  wavelength/Rossby  radius  ratio,  amplitude/wavelength  ratio,  and 
cross  stream  Froude  number.  Tests  were  run  using  oceanographic  con¬ 
ditions  representative  of  two  regions,  the  Gulf  Stream  downstream  of 


Cape  Hatteras,  and  the  Sargasso  Sea  at 


30^5. 


The  model  was  shown  to  be  sensitive  to  the  wavelength/Rossby 
radius  ratio.  For  the  Sargasso  Sea  conditions  a  30  km  disturbance  was 
more  stable  with  respect  to  amplitude  growth  than  a  100  km  disturbance. 
Both  the  magnitude  and  sign  of  the  cross  stream  Froude  number  affected 
the  meander  stability,  with  a  negative  Froude  number  (downward  entrain¬ 
ment  near  the  surface  front)  showing  the  greatest  amplitude  growth. 

The  model  results  were  very  sensitive  to  the  meander  amplitude.  Ampli¬ 
tude  growth  rates  obtained  from  this  model  for  small  amplitude  meanders 


■were  similar  to  results  from  several  small  amplitude  analytic  models. 
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The  growth  rates  from  large  amplitude  tests  were  much  lower,  and 
closely  matched  observational  results. 

This  model  also  showed  a  tendency  toward  growing  inertial 
period  oscillations  near  the  boundary  between  the  dissipative  and 
inviscid  regions.  These  oscillations  were  controlled  here  by  inter¬ 
facial  friction,/  A  more  accurate  representation  of  the  dissipative 

I  \ 

region  may  be  necessary  to  properly  dissipate  the  inertial  period 

energy.  \ 
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INTRODUCTION 


1.1  -  Observations  of  Oceanic  Fronts 

Oceanic  fronts  constitute  a  boundary  between  two  water  masses  of 
different  temperature,  salinity,  and  density.  The  change  between  the 
water  masses  may  occur  across  a  distance  of  a  few  meters  to  a  few  kilo¬ 
meters.  The  front  generally  angles  downward  with  a  shallow  slope  (e.g., 
one  in  one  hundred).  As  suggested  by  Stommel  (1976),  this  downward 
slope  often  tapers  off  to  give  an  'exponential'  profile.  The  abrupt 
change  in  physical  characteristics  across  the  front  has  been  observed  to 
have  a  profound  effect  on  the  dynamics  of  the  region  (Voorhis,  1969; 
Stommel,  1976),  and  upon  the  biota.  Backus  et  al.  (1969)  found  major 
changes  in  the  biota  across  a  front  in  the  Sargasso  Sea.  The  species 
composition,  total  fish  population,  and  primary  productivity  were  all 
strikingly  different  on  the  two  sides  of  the  boundary.  A  cross-frontal 
shear  in  the  velocity  component  parallel  to  the  front  is  often  observed, 
and  is  often  quite  intense.  In  addition,  significant  surface  velocities 
perpendicular  to  the  front  have  been  observed,  and  the  front  is  a  zone 
of  surface  convergence  and  sinking  (Garvine  and  Monk,  1974;  Voorhis, 
1969).  This  convergence  collects  nutrients,  debris,  and  foam  at  the 
surface  expression  of  the  frontal  interface,  often  making  it  highly 
visible.  Other  characteristics  of  the  water  such  as  color,  turbidity, 
and  surface  roughness  may  change  abruptly,  thus  also  enhancing  the 
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visibility  of  the  boundary. 


Studies  of  fronts  in  the  ocean  and  atmosphere  can  be  divided 
into  three  groups,  depending  upon  the  stage  in  the  evolution  of  the 
front  which  is  examined.  There  are  those  studies  which  are  concerned 
with  the  generation  of  fronts,  or  frontogenesis;  those  which  deal  with 
the  behavior  of  established,  persistent  fronts;  and  those  involving  the 
dissipation  and  disappearance  of  fronts,  or  frontolysis.  For  many  types 
of  fronts,  the  frontogenic  and  frontolytic  periods  occur  in  such  close 
apposition  that  a  middle  period  of  enduring  existence  cannot  be  defined. 
However,  other  types  of  fronts,  such  as  the  Gulf  Stream  and  the  atmo¬ 
spheric  jet  stream,  have  such  a  degree  of  permanence  that  it  is  the 
generative  and  dissipative  periods  which  do  not  exist.  These  semi¬ 
permanent  fronts  can  have  very  active  dynamics  and  often  strongly 
influence  the  local  circulations  and  energy  exchanges. 

A  major  feature  of  persistent  fronts  is  their  tendency  to 
develop  meanders,  i.e.,  displacements  of  the  frontal  interface  away  from 
a  straight  line.  Meanders  are  most  often  observed  as  lateral  displace¬ 
ments  in  the  position  of  the  surface  expression  of  the  front,  forming  a 
wavelike  pattern.  In  western  boundary  currents  such  as  the  Gulf  Stream, 
meanders  can  grow  to  such  a  large  amplitude  that  the  bulge  breaks  off  to 
form  a  cyclonic  or  anticyclonic  eddy  or  ’ring’.  South  of  Cape  Hatteras 
the  Gulf  Stream  usually  exhibits  relatively  small  amplitude  fluctuations. 
After  the  Stream's  departure  from  the  shelf/slope  region  north  of  Cape 
Hatteras,  however,  the  meanders  increase  in  amplitude  and  often  form 
eddies  south  of  Georges  Bank. 


3 


Some  of  the  earliest  observations  of  Gulf  Stream  meanders  were 
taken  during  the  Multiple  Ship  Survey  of  June  1950  (Fuglister  and 
Worthington,  1951).  More  recent  shipborne  bathythermograph  data  from 
Operation  Cabot  (Hansen,  1970)  resulted  in  estimates  for  the  dominant 
wavelength  and  phase  speed  of  320  km  and  8  cm/sec.  These  results  were 
confirmed  by  Halliwell  and  Mooers  (1978)  in  a  study  using  satellite 
thermal  imagery.  This  study  found  the  dominant  meander  mode  to  have  a 
wavelength  of  320  km,  a  phase  speed  of  6  cm/sec,  and  a  frequency  of  6 
cycles  per  year.  Meanders  have  also  been  observed  along  fronts  in  the 
Sargasso  Sea  by  Voorhis  (1969)  and  Katz  (1969)  and  along  the  front  of 
the  California  Current  by  Bernstein  et  al.  (1977). 

1.2  -  Frontal  Models 

Several  theoretical  models  have  been  developed  to  study  the 
motions  of  established  fronts,  particularly  for  the  atmospheric  case. 
Modeling  of  oceanic  fronts  has  not  progressed  as  rapidly,  due  in  part  to 
the  difficulty  in  verification  and  to  the  more  complex  environs  of  the 
ocean.  An  important  difference  between  the  atmosphere  and  the  ocean  is 
that  the  atmospheric  density  is  controlled  almost  entirely  by  temperature, 
while  the  oceanic  density  is  controlled  by  both  temperature  and  salinity. 
However,  several  atmospheric  models  capture  much  of  the  essential  physics 
of  oceanic  fronts,  so  that  they  are  relevant  to  this  discussion. 

Stoker  (1953)  rigorously  formulated  the  full  nonlinear  problem 
for  the  propagation  of  an  atmospheric  frontal  disturbance.  He  set  up 
four  models  with  four  levels  of  complexity,  from  few  restrictions  to 


highly  restricted.  However,  even  the  least  complex  (most  restricted) 
model  was  inherently  nonlinear  and  none  of  the  models  could  be  solved 
analytically.  Since  general  purpose  computers  were  not  available  at 
that  time,  a  numerical  approach  was  not  practical  either,  so  no  solution 
was  obtained.  The  assumptions  made  for  the  most  restricted  model 
included  an  inviscid,  incompressible  fluid,  no  change  in  the  Coriolis 
parameter  with  latitude  (f-plane  approximation),  initially  geostrophic 
balance  across  the  frontal  interface  which  was  assumed  to  be  a  planar 
discontinuity,  a  continuous  pressure  field,  vertical  hydrostatic  balance, 
vanishing  fluid  velocity  relative  to  the  front  perpendicular  to  the 
front  both  at  the  front  and  well  away  from  the  front,  and  a  cold  (lower) 
layer  much  thinner  than  the  ambient  fluid. 

Whitham  (1953),  in  a  companion  paper  approached  Stoker’s  most 
restricted  model  using  the  method  of  characteristics,  which  permits  a 
graphic  type  of  solution.  His  solution  showed  the  nonlinear  character 
of  disturbance  propagation  and  the  tendency  toward  cyclogenesis. 

Kasahara  et  al.  (1965)  solved  Stoker's  problem  III  (slightly  less 
restricted  than  the  problem  solved  by  Whitham)  numerically.  The  problem 
was  solved  as  a  free  boundary,  initial  value  problem  by  using  a  finite 
difference  scheme.  An  initially  sinusoidal  disturbance  was  shown  to 
propagate  along  the  front  from  west  to  east  and  to  result  in  incipient 
cyclogenesis. 

To  achieve  solutions  to  the  full  nonlinear  problem,  both  Whitham 
and  Kasahara  et  al.  used  equations  which  were  a  rough  approximation  to 
the  physical  system  and  then  +  -ied  to  achieve  exact  solutions.  If, 


however,  an  approximate  solution  to  the  equations  is  sought,  the  equa¬ 
tions  describing  the  physical  system  may  be  placed  under  lighter  restric 
tions.  The  most  common  approximation  technique  is  perturbation  analysis 
This  method  has  been  very  successful  in  many  fields  of  theoretical  and 
applied  physics.  The  assumption  which  is  required  for  the  use  of  this 
method  is  that  the  disturbance  is  a  small  perturbation  on  the  mean  flow 
(amplitude  much  smaller  than  the  wavelength).  This  allows  the  solution 
to  be  expanded  in  a  series,  from  which  the  second  order  and  higher  terms 
are  dropped,  thus  leaving  a  linear  problem,  rather  than  a  nonlinear  one. 
Often  used  in  conjunction  with  perturbation  methods  is  Fourier  analysis 
in  which  the  disturbance  is  assumed  to  be  wavelike  and  sinusoidal.  This 
method  is  of  particular  utility  when  the  frequencies  cannot  be  chosen  a 
priori. 

Models  which  rely  on  the  perturbation  method  of  solution  must 
assume  that  the  meander  amplitudes  are  small.  However,  prior  to  occlu¬ 
sion,  the  amplitudes  may  exceed  the  wavelength,  and  in  any  case  are 
certainly  not  small.  Perturbation  models  exclude  the  nonlinear  inter¬ 
actions  which  must  be  occurring  in  the  frontal  regions  by  dropping 
higher  order  terms  in  the  series  expansion.  While  small  amplitude 
models  may  well  represent  the  early  stages  of  meander  propagation,  they 
necessarily  omit  much  of  the  essential  physics. 

One  of  the  first  studies  to  use  perturbation  analysis  for  the 
study  of  oceanic  fronts  was  that  of  Ouxbury  (1963).  Duxbury's  system 
included  a  planar  front  with  a  constant  velocity  shear  across  it.  The 
model  was  inviscid  with  no  cross-stream  motion  in  the  basic  state. 
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All  solutions  were  dynamically  stable,  due  largely  to  the  neglect  of  the 
coupling  with  the  dynamics  in  the  lower  layer,  as  shown  by  Qrlanski 
(1969).  The  degree  of  stability  was  shown  to  be  dependent  upon  the 
Froude  number.  Orlanski  (1968,  1969)  derived  perturbation  models  that 
did  show  unstable  wave  modes.  The  flow  of  energy  in  these  models  was 
shown  to  be  from  the  mean  potential  to  the  mean  kinetic  energy  compo¬ 
nent,  as  expected  from  large  scale  turbulence  and  energy  considerations. 
However,  neither  of  Orlanski 's  models  included  the  effects  of  dissi¬ 
pation  or  cross-stream  flow  in  the  initial  state. 

Hansen  (1970)  made  observations  of  Gulf  Stream  meanders  by 
shipborne  bathythermograph  and  compared  his  observational  results  with 
the  model  results  of  several  authors  (Ouxbury,  1963;  Haurwitz  and 
Panofsky,  1950;  Lipps,  1963;  Tareev,  1965;  and  Orlanski,  1969).  The 
results  showed  that  the  barotropic  models  (Haurwitz  and  Panofsky,  Lipps, 
and  Duxbury)  were  not  accurate  in  selecting  the  wavelength  of  the 
dominant  mode.  Only  the  Orlanski  model  made  a  good  estimate  of  the 
phase  speed  and  only  the  models  of  Lipps  and  Tareev  estimated  the  rate 
of  meander  spatial  growth  well.  None  of  the  estimates  of  temporal 
growth  were  particularly  accurate.  Thacker  (1976)  showed  that,  de¬ 
pending  upon  the  mean  velocity  and  the  velocity  shear  across  the  inter¬ 
face,  the  Gulf  Stream  can  show  either  predominantly  spatial  or  predomi¬ 
nantly  temporal  growth. 

Since  these  models  made  consistently  poor  estimates  of  certain 
important  quantities,  some  fundamental  phys<cs  must  be  missing  in  their 
formulations.  Possibilities  for  improvement  include:  the  use  of  a 
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dissipation  mechanism  such  as  friction,  cross-stream  flow  in  the  basic 
state,  and  the  easing  of  the  small  amplitude  restriction. 

In  studies  of  large  scale  oceanic  processes  the  fluid  is  often 
assumed  to  be  inviscid.  However,  this  may  not  be  valid  for  frontal 
dynamics.  Grammeltvedt  (1970)  in  a  numerical  solution  of  the  problem 
considered  by  Kasahara  et  al.  (1965)  (and  proposed  by  Stoker,  1953) 
showed  the  development  of  asymmetries  leading  to  occlusion  and  eddy 
formation.  However,  to  maintain  numerical  stability,  the  artifice  of  a 
dissipative  mechanism  had  to  be  introduced.  Assaf  (1977),  in  an  exten¬ 
sion  of  an  earlier  model  by  Charney  for  the  Florida  Current,  produced  a 
cyclonic  shear  zone  by  using  vertical  (interfacial)  friction  rather  than 
lateral  friction.  He  also  showed  that  friction  may  act  to  increase  the 
available  potential  energy  in  the  Sargasso  Sea.  It  is  also  apparent 
that  friction  generated  by  a  velocity  shear  at  a  boundary  would  be  more 
important  where  the  fluid  layer  is  thin  than  where  it  is  deep.  Thus,  it 
appears  that  some  dissipative  mechanism  (such  as  friction)  is  important 
to  the  dynamics  of  the  front,  particularly  near  the  surface  front. 

To  show  the  importance  of  cross-frontal  motions,  Saltzman  and 
Tang  (1975)  derived  an  analytical  model  of  an  atmospheric  front.  They 
showed  that  in  the  upper  troposphere  (mid-depth  in  an  oceanic  front) 
there  is  a  horizontal  convergence  upstream  and  a  divergence  downstream 
of  a  meander  trough,  resulting  in  a  non-zero  cross  front  flow.  The 
results  indicate  that  second  order,  ageostrophic  effects  can  modify  a 
two  layer  flow.  The  implication  is  that  asymmetric  features  such  as 


frontal  meanders  are  the  result  of  a  large  scale  baroclinic  instability 
with  second  order  effects  superposed  near  the  surface  front.  Thus, 
there  may  be  two  flow  regimes,  one  encompassing  the  entire  frontal 
region  and  one  active  only  near  the  surface  front.  In  another  atmo¬ 
spheric  frontal  model,  Orlanski  and  Ross  (1977)  showed  that  a  cross¬ 
stream  circulation  could  be  caused  by  an  ageostrophic  residue  from  the 
vertical  shear  of  the  horizontal  velocity  and  the  horizontal  gradient  of 
the  temperature.  In  this  model  the  cross-stream  motion  was  required  so 
that  the  steady  state  would  be  maintained.  Newton  (1978)  showed  that 
mid- tropospheric  fronts  (such  as  the  jet  stream)  have  more  similarities 
to  oceanic  fronts  than  do  the  lower  tropospheric  fronts  which  intersect 
the  ground.  He  discussed  analogies  between  mid-tropospheric  and  oceanic 
fronts  and  derived  a  model  for  the  effects  of  confluence  on  the  circu¬ 
lation  of  an  oceanic  frontal  region.  The  concentration  of  horizontal 
gradients  due  to  confluence  in  meander  motions  was  shown  to  cause  a 
sinking  at  mid-depths  and  a  convergence  at  the  surface  in  meander  troughs, 
so  that  the  isopycnals  become  closer  in  the  troughs.  This  agrees  with 
the  results  of  Saltzman  and  Tang  (1975)  and  is  supported  by  limited 
observation  for  atmospheric  fronts.  Formato  (1979),  in  a  model  which 
explicitly  includes  frictional  dissipation  near  the  surface  front, 
showed  that  the  cross-stream  flow  is  non-zero  both  in  the  steady  state 
and  in  the  presence  of  meanders.  Thus,  cross- stream  motions  are  an 
imoortant  and  perhaps  essential  component  of  the  frontal  dynamics. 

Garvine  (1974)  derived  a  model  for  a  small  scale,  steady  state 
oceanic  front  which  included  friction  and  mass  entrainment.  These 
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mechanisms  were  shown  to  be  required  at  small  scales  for  the  maintenance 
of  a  steady  state.  This  model  was  extended  (Garvine  1979a;  1979b)  to 
include  the  rotational  (Coriolis)  effects  present  in  larger  scale  fronts. 
This  model  was  hydrodynamic  (density  field  not  determined  by  model 
orocesses),  vertically  integrated,  and  used  bulk  parameterizations  for 
quantities  such  as  friction.  An  important  result  of  these  papers  was 
the  derivation  of  the  rotation  parameter,  which  expresses  the  relative 
importance  of  the  rotational  effects  to  the  effects  of  friction.  Small 
values  indicate  that  frictional  effects  are  dominant.  A  value  of  zero 
reduced  the  1979  model  to  the  1974  case.  For  large  scale  fronts,  such 
as  the  Gulf  Stream,  the  rotation  parameter  is  large  so  that  rotational 
effects  are  dominant  for  the  front  as  a  whole.  Garvine  demonstrated 
that  dissipation  processes  are  important  only  for  a  distance  of  10%  to 
20%  of  the  baroclinic  Rossby  radius  (the  cross-stream  scale  of  the 
frontal  system)  away  from  the  surface  front.  However,  due  to  mass 
entrainment  and  mixing  in  this  area  the  entire  frontal  zone  is  affected 
by  a  cross-stream  mass  flux. 

To  improve  upon  present  models,  three  factors  should  be  con¬ 
sidered:  a  dissipative  mechanism  such  as  friction  which  is  active  near 
the  surface  expression  of  the  front,  motion  perpendicular  to  the  front 
in  the  upper  (lighter)  layer,  and  the  nonlinear  effects  of  finite 
amplitude  meanders.  The  Intent  of  this  paper  is  to  present  a  model  of 
an  upper  ocean  density  front  which  specifically  addresses  the  three 
aforementioned  considerations.  This  model  is  based  on  the  work  of 
Kasahara  et  al.,  but  is  altered  according  to  the  results  of  Garvine 
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and  Formato.  The  Kasahara  et  al.  atmospheric  front  model  allows  for  the 
full  nonlinear  equations,  and  the  modifications  introduced  here  add 
dissipation  near  the  surface  front  and  cross  stream  mass  flux. 

The  present  model  is  a  conjunction  of  two  models.  A  fixed  grid 
Eulerian  model  is  used  for  the  main  portion  of  the  flow  and  is  con¬ 
strained  by  three  fixed  boundaries  at  which  standard  boundary  conditions 
are  applied.  However,  the  fourth  boundary,  the  front,  must  be  free,  a 
condition  which  an  Eulerian  model  cannot  deal  with  directly.  To  accurately 
compute  the  behavior  of  the  front,  a  second,  Lagrangian,  model  is  used 
along  the  front.  The  Eulerian  grid  space  provides  the  acceleration 
terms  for  the  Lagrangian  model  through  the  gradient  of  the  interfacial 
depth.  The  results  along  the  free  boundary  are  used  at  each  time  step 
to  form  a  boundary  condition  for  the  Eulerian  model.  The  two  models 
thus  complement  each  other  and  are  integrated  in  parallel  through  time 
to  describe  the  time  evolution  of  the  flow. 


1 


THE  EULERIAN  MODEL 


2.1  -  Physical  Description 

To  model  an  upper  ocean  density  front  we  assume  an  incompres¬ 
sible,  hydrostatic  fluid  on  a  rotating  earth.  The  flow  is  assumed  to 
be  periodic  along  the  direction  of  the  primary  flow.  The  periodicity 
is  imposed  by  the  boundary  conditions  and  the  length  scale  is  thus  equal 
to  the  length  of  the  model  grid.  The  upper  surface  of  the  fluid  is  flat 
and  rigid.  A  two  layer  system  is  used  with  the  upper  layer  divided  into 
two  regions,  dissipative  and  inviscid.  Figure  1  shows  a  vertical 
section  across  the  front  to  illustrate  the  division  of  the  flow  field. 
The  lower  layer,  with  the  greatest  volume,  is  called  the  ambient  water 
and  lies  below  and  to  the  left  of  the  frontal  (pycnocline)  interface. 

The  characteristic  depth  is  h  and  the  density  is  taken  to  be  uniform  at 
o,,.  The  other  two  regions  are  comprised  of  lighter  water  and  lie  above 
the  interface.  One  region  is  the  outer  (inviscid)  zone  where  the 
density  is  uniform  at  pm  -  Ap.  The  depth  of  the  interface  is  D,  and 
varies  in  space  and  time.  Since  the  flow  in  this  region  has  been 
taken  to  be  inviscid,  there  is  no  vertical  shear  in  the  horizontal 
velocity  components,  u  and  v.  To  the  right  of  the  inviscid  region  of 
the  front  is  a  reservoir  of  light  water,  called  the  parent  pool,  where 
the  depth  and  velocity  remain  unaffected  by  any  variations  in  the 
motion  in  the  frontal  zone.  The  horizontal  scale  of  the  entire  frontal 
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zone  is  the  baroclinic  Rossby  radius,  \  -  c/f.  Here  c  is  the  internal 
wave  speed,  c  =  (g'D0)1/2,  where  g*  =  g,  the  reduced  gravity  and  DQ 

so 

is  the  depth  of  the  pycnocline  in  the  parent  pool.  The  Coriolis  parameter 
is  f,  and  is  taken  to  be  constant  with  latitude  (f-plane  approximation). 

The  second  region  above  the  interface  is  the  dissipative  zone. 

This  region  adjoins  the  surface  front  and  extends  away  from  it  for  a 
distance  of  10%  to  20%  of  \.  The  primary  importance  of  the  inner  region 
is  its  influence  on  the  outer  region  at  the  common  juncture.  The 
assumed  high  level  of  turbulence  in  the  inner  region  permits  the  use  of 
a  local,  quasi-steady  solution,  because  the  time  scale  there  for  flow 
adjustments  to  changes  in  the  boundary  conditions  is  short  compared  to 
the  time  scale  for  the  outer  region.  In  the  model  the  inner  (dissi¬ 
pative)  region  is  treated  as  a  'grey  box'.  Some  of  the  details  of  the 
turbulent  processes  occurring  inside  of  it  are  known,  but  not  enough  for 
it  to  be  modeled  accurately  at  a  level  compatible  with  this  modeling 
effort.  Therefore,  the  inner  region  is  assumed  to  be  a  source  or  sink 
of  water  for  the  outer  region,  but  the  details  of  the  flow  are  not 
modeled. 


We  assume  that  D/h  is  small  which  permits  great  simplification 
in  the  dynamics,  since  the  motion  in  the  ambient  layer  is  then  uncoupled 
from  that  above  the  interface  (Garvine,  1979a),  and  will  be  assumed  to 
be  in  uniform  geostrophic  balance.  This  means  that  above  the  pycnocline, 
including  the  dissipative  zone,  the  flow  will  be  in  isostatic  balance. 
This  also  allows  us  to  reduce  the  computation  to  a  one  layer,  reduced 
gravity  model. 
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Figure  2  shows  the  coordinate  system  to  be  used  in  the  hori¬ 
zontal  plane.  In  the  basic  or  unperturbed  state,  the  surface  front  lies 
parallel  to  the  y-axis.  The  positive  x-axis  points  towards  the  light 
water  reservoir,  and  the  z-axis  is  upwards.  This  coordinate  system  is 
similar  to  that  used  by  Formato.  Here,  however,  the  perturbations  to 
the  front  will  not  be  decomposed  into  steady  state  and  perturbation 
components.  A  horizontal  perturbation  may  be  applied  to  the  inner/outer 
boundary  with  an  amplitude  to  wavelength  ratio  between  zero  (straight  or 
unperturbed)  to  order  one.  The  perturbation  is  applied  as  a  cosine 
shaped  disturbance. 

2.2  -  Initial  Equations 

The  system  of  equations  governing  the  inviscid  region  are  well 
known.  Listed  i  order  of  vertically  integrated  continuity,  x  momentum, 
and  y  momentum  they  are,  for  an  earth  fixed  coordinate  system: 


0t*  +  (u*0)x*  +  (v*D)y* 

8  0 

(2.2.1) 

'$*  +  +  V*uj*  -  fv*  8 

p*  px* 

(2.2.2) 

t*  +  u*v£*  +  V*VJ*  +  V  3 

p  py* 

(2.2.3) 

These  equations  are  the  same  as  those  for  large  amplitude 
(nonlinear)  long  internal  waves,  including  the  effects  of  the  earth's 
rotation.  For  the  non- rotational  case  (f  8  0),  they  are  identical  to 
the  nonlinear  shallow  water  wave  equations,  with  0  analogous  to  the 
surface  height.  The  system  is  thus  hyperbolic  and  has  real  character 
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i sties.  As  is  well  known,  such  wave  systems  are  often  found  in  nature 
with  hydraulic  jumps  at  their  advancing  boundaries  (Stoker,  1957). 

Tidal  bores  are  a  common  example.  In  this  case  the  dissipative  region 
is  directly  analogous  to  the  shallow  water  bore  and  the  inviscid  region 
is  analogous  to  the  long  wave  region  behind  the  bore. 

For  generality  we  include  a  transformation  to  a  coordinate 
system  which  is  moving  with  the  speed  of  the  unperturbed  portion  of  the 
front.  This  will  allow  the  separation  of  the  linear  frontal  propagation 
from  the  frontal  meandering.  A  simple  Galilean  transformation  is 


x*  =  x*  -  Uft*  ’ 


y'  =  y*  -  vft* 


t'  =  t* 


u1  =  u*  -  U, 


v 1  *  v*  -  V, 


(2.2.4) 


where  and  V*  are  the  frontal  propagation  velocities.  This  leads  to 


0t,  +  (u'D)x.  +  (v'D)  ,  =  0 


ut'  +  u'ux'  +  v'uy'  "  f(v'  +  V  =  "  o  Px' 


ui,  +  u’v'  +  v'v\  +  f(u'  +  U£)  = 


'  o  py' 


(2.2.5) 


(2.2.6) 


(2.2.7) 
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The  pressure  terras  on  the  right  hand  side  will  be  difficult  to 
work  with  in  this  type  of  model.  An  easier  variable  to  use  is  the 
interfacial  depth.  Since  the  upper  layer  is  homogeneous  and  in  isostatic 
balance,  this  conversion  may  be  made  by  calculating  the  pressure  field 
hydrostatically. 


P  *  g(p„  -  Ap)(n  -  z) 

-D  <  z  <  o 

(2.2.3) 

P  =  ^(p,,  -  Ap)n  -  (p^z  +  ApD)} 

2  <  -D 

(2.2.9) 

Neglecting  terras  of  order  compared  to  unity,  the  hydrostatic 

pressure  gradient  is 

7h  t  3  «V 

00 

-D  <  2  <  n 

(2.2.10) 

7h  !r  *  *V("  -  f- 

oo  ao 

2  <  -D 

(2.2.11) 

Below  the  interface,  in  the  ambient  fluid  region,  the  hydro¬ 
static  pressure  gradient  is  balanced  geostrophically.  Setting  U, ,  V. 

a  a 

as  the  ambient  fluid  velocities,  we  have 


g,h(„  -  f-  0)  .  f  va  t  -  f  ua  T 

CO 

Performing  the  Galilean  transformations  as  above 


(2.2.12) 


(2.2.13) 


1 


so  that 

gvhn  -  9  ^  VhD  =  f(vM  +  Vf)t  -  ftu^  -  Uf)j 

Using  the  reduced  gravity  g'  =  g  above  the  interface 

%  f-  -  g^n  •  g'V  t  f(v.  *  vf)T 
00 

*  +  Uf)T 

so  that  (2.2.6)  and  (2.2.7)  become 

ut'  +  u'ux*  +  v<  uy 1  ■  f(v‘  +  vf)  =  -9°x' 

-  fv  -  fV£ 

m  T 


v|,  +  u'v',  +  v'v',  +  f(u'  +  Uf)  =  -gO  , 
"  *  y  T  y 


-  fu  +  fu. 


The  basic  equations  for  the  system  are 


(2.2.14) 


(2.2.15) 


(2.2.16) 


(2.2.17) 


+  u'u',  +  v'u',  -  fv'  =  -g‘Dv,  -  fv' 
a  y  x  » 

(2.2.18) 

+  u'v',  +  v'v',  +  fu'  =  -g'D  .  +  fu' 
x  y  y  * 

(2.2.19) 

°t'  +  (u'0)x'  +  (v'D)y'  =  0 

(2.2.20) 
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2.3  -  Scaling 

For  convenience  the  basic  equations  (2.2.18)-(2.2.20)  will  now 
be  put  into  non-dimensional  form.  The  primary  assumptions  made  in  this 
scaling  are  that  the  velocity  scale  is  the  internal  wave  speed  C  -  (g ' DQ) 
and  the  length  scale  is  the  baroclinic  Rossby  radius  of  deformation 
x  *  j.  The  obvious  time  scale  is  then  based  on  the  inertial  frequency, 
f.  With  these  scales,  the  variable  transformations  are 


X  = 

Xn 

u* 

=  UAf 

y  a 

Xz 

v' 

»  vxf 

t  = 

t/f 

u* 

=  u  xf 
00 

o 

it 

5Do 

v' 

CO 

3  V^Xf 

g'D  *  $(fx)2 


Note  that  the  internal  cross  stream  Froude  number  F„  *  -~ 

r  c 

which  is  equal  to  uQ>  the  parent  pool  cross  stream  velocity. 

Using  these  equations  (2.2.18)  becomes 


f  ~  (ufx)  +  ufx  (ufx)  +  vfx  (ufx) 


-  f(vfx)  =  -  (3f2X2)  -  f(vjx)  (2.3.2) 

Multiplying  by  -i-  yields 

rx 
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(2.3.3) 


U  +  UU  +  VU  -  V  =  -5  -  V 

t  n  ?  n 

Equation  (2.2.19)  is,  similarly 

v  +  uv  +  vv  +  u  3  -<$  +  u  (2.3.4) 

t  n  ?  n  « 

With  the  new  variables  equation  (2.2.20)  is 

5  =  -  u6  -  Su  -  v<5  -  5v  (2.3.5) 

t  n  n  t  z 

2.4  -  Conservation  form  of  the  equations 

The  difference  scheme  to  be  used  requires  that  the  initial 
equations  be  in  conservation  form  (see  appendix  I).  This  means  that 
they  must  be  in  the  form 


wt  *  fx  +  9y  +  z  (2.4.1) 

where  w  is  a  vector  containing  the  three  dependent  variables.  The 
vectors  f  and  g  are  functions  of  w  which  are  differentiated  in  space, 
and  z  is  a  vector  containing  the  inhomogeneous  terms.  Equation  (2.3.5) 
is  already  in  this  form  and  equations  (2.3.3)  and  (2.3.4)  may  be  trans¬ 
formed  by  converting  to  mass  flux  variables  instead  of  the  velocity 
terms.  Multiplying  (2.3.3)  by  S 

Su  +  <5uu  +  <5vv  -  SV  3  - 66  -  5v  (2.4.2) 

t  n  s  n  * 

and  multiplying  (2.3.3)  by  u 

2 

u<5  +  u  5  +  <$uu  +  uv6,  +  Suv„  =  0 

x  n  n  C  C 


(2.4.3) 


and  adding  the  two 


6u  +  U6  +  6UU  +  Suu  +  u  0  +  66  +  UV6„ 

t  t  n  n  n  n  5 


+  6UV„  +  5VU  -  iSV  +  6V  =0 
5  5  00 


w,  *  -M  -[?*%- 

n  V  1 


+  6v  -  6v  (2.4.4) 


Similarly  we  may  multiply  (2.3.4)  by  6  and  (2.3.5)  by  v  and  add 


the  two  to  give 


oV  +  V5  +  (SUV  +  (SVU  +  UV6  +  25VV 
t  r  n  n  n  n 


+  V  <5,  +  66  +  6u  -  oU  *  0 

5  5 


<**),  •-(¥*)  -  M 


-  6u  -  5u  (2.4.5) 


If  we  let  <t>  *  5u,  and  t|>  =  sv  equations  (2.4.4),  (2.4.5),  and 


(2.3.5)  become 


=•  -f4-+  1  621  +  +  *  - 


<p  -  5v 


(2.4.6) 
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=  - 


2  1 

U  j 

n 

6  2  0 
,  . 

-  $  +  6U 


(2.4.7) 


(5 

T 


(2.4.8) 


or 


where 


wT  s  f(w)n  +  g(w);  +  2 


(2.4.9) 


> 

i  ! 

: 

f  1_  .  1  K  i 
6  +?5 

w  =  < 

f  =  - 

6 

.  « 

* 

6 

-  SV 

00 

g  = 

ii  +  i  .2 

■  z  =  • 

-4>  +  <$u 

T  00 

.  0 

(2.4.10) 


This  is  the  form  of  the  equations  that  will  be  inserted  into  the 
finite  differencing  scheme  for  the  fixed  grid  portion  of  the  flow. 


2.5  -  Interfacial  friction 


Although  the  outer  region  is  modeled  as  an  inviscid  fluid  there 
will  be  some  frictional  effects,  primarily  at  the  interface  between  the 
water  in  the  outer  region  and  the  ambient  water.  This  friction  may  be 
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characterized  by  a  bulk  stress  formula 

Fq  =  Cq  V | v* |  (2.5.1) 

where  Fg  is  the  frictional  drag  force  and  Cg  is  the  coefficient  of 
friction.  Note  that  Fg  is  proportional  to  the  square  of  the  velocity, 
although  the  absolute  value  form  is  used  here  to  orient  the  drag  force 
properly  with  respect  to  the  velocity.  The  coefficient  of  drag  must  be 
scaled  in  the  same  manner  as  the  main  equations  were  in  section  2.3. 
Using  the  velocity  scaling  we  have: 

Foxf2  55  cd  r  vMx2f2 


or 


Fg  =  Cg  v | v |  (2.5.2) 

The  use  of  interfacial  drag  provides  a  simple  method  for  in¬ 
cluding  the  effect  of  local  instabilities  of  high  frequency  internal 
gravity  waves  propagating  along  the  interface.  It  also  stabilizes  the 
model  by  reducing  the  inertial  oscillations  which  occur  as  a  result  of 
imperfect  initial  conditions  and  which  can  amplify  if  not  suppressed. 

The  friction  in  this  model  is  not  treated  as  operating  within 
the  upper  layer,  but  only  at  the  interface,  and  is  described  by  a  bulk 
parametrization.  Therefore,  it  is  not  included  in  the  initial  momentum 
equations  (2. 2. 2-2. 2. 3)  and  is  thus  excluded  from  the  results  of  the 
finite  difference  computations.  The  calculation  of  friction  is  appended 
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to  the  differencing  as  an  auxiliary  operation,  along  with  the  adjust¬ 
ments  for  the  boundary  conditions  discussed  in  sections  2.3  and  2.9. 

2.6  -  Kinetic  and  potential  energy 

The  behavior  of  the  model  in  the  outer  region  will  be  followed 
by  the  time  and  space  evolution  of  the  three  independent  variables; 
depth  and  the  two  components  of  velocity.  The  spatially  integrated 
kinetic  and  potential  energies  will  also  be  computed  to  assist  in  the 
physical  interpretation  of  the  model  results.  The  equations  for  the 
kinetic  and  potential  energies  spatially  sunuied  over  the  outer  region 
are: 


K.E.  =  \  H  (u*2  +  vl2)D  (2.6.1) 

n  C 

P.E.  =  £g'  H  D2  (2-6.2) 

n  c 

Scaling  these  as  in  section  2.3: 

K.E.  *  \  D0  c2  H  (u2  ♦  v2)5  (2.6.3) 

n  c 

P.E.  *  ^f)2  l  l  52  (2.6.4) 

n  c 

Following  Qrlanski  and  Cox  (1973)  these  are  divided  into  mean  and  per¬ 
turbation  components.  Horizontal  means  are  taken  in  the  along  frontal 
direction  for  each  of  the  three  basic  variables  as: 
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5(x,t)  =  7T  y(x,y,t)dy 


(2.6.5) 


and  the  perturbation  quantities  as: 


u'(x,y,t)  =  u(x,y,t)  -  u(x,t) 


(2.6.6) 


From  the  mean  and  perturbation  components  of  the  three  independent 
variables,  the  mean  and  perturbation  components  of  the  kinetic  and 
potential  energies  are  computed  by  equations  (2.6.3)  and  (2.6.4).  The 
results  of  these  calculations  will  indicate  the  directions  and  relative 
rates  of  the  energy  transformations.  These  will  be  diagnostic  of  the 
stability  characteristics  of  the  perturbations  to  the  flow. 


2.7  -  Potential  vorticity 


In  the  investigation  of  wave  phenomena  it  is  sometimes  useful  to 
employ  a  vorticity  equation.  Such  an  equation  may  be  easily  derived  for 
this  case.  By  cross  differentiating  (2.2.18)  and  (2.2.19)  and  then 
subtracting,  the  terms  involving  5  may  be  eliminated.  Setting  the 
vorticity  as  y  =  v^  -  u^,  we  have 


Yt  +  (UYn)  +  (VY)?  +  un  +  V  ■  0 


(2.7.1) 


Y  +  (uy  +  U)  +  (VY  +  V)  ■  0 
t  n  c 


(2.7.2) 


By  combining  this  equation  with  continuity  (2.2.20),  we  may 
derive  an  equation  for  the  conservation  of  potential  vorticity: 
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(2.7.3) 


In  analytical  models  such  an  equation  is  often  substituted  for 
one  of  the  momentum  equations  in  the  set  of  initial  equations.  In  this 
case  the  vorticity  equation  does  not  fit  into  the  differencing  scheme 
and  attempts  to  circumvent  the  problem  by  a  spatial  integration  were 
unsuccessful . 

The  potential  vorticity  equation  is  still  useful,  however,  in 
checking  the  results  of  the  model.  The  potential  vorticity  at  each  grid 
point  may  be  calculated  from  (2.7.2)  and  displayed  graphically.  Since 
the  potential  vorticity  field  is  originally  uniform  from  the  initial 
conditions,  and  since  the  outer  region  is  inviscid,  except  as  discussed 
in  section  2.5,  the  uniformity  of  the  field  at  later  points  in  time  is  a 
test  of  the  stability  and  accuracy  of  the  model. 

2.3  -  Conditions  on  the  end  boundaries 

Since  the  finite  difference  scheme  can  cover  only  a  limited 
region  of  space  along  the  front,  boundary  conditions  must  be  applied 
along  lines  perpendicular  to  the  front  at  the  ends  of  the  grid  space 
(the  'top'  and  'bottom'  boundaries  in  Fig.  3).  It  is  assumed  that  these 
edges  are  far  away  relative  to  the  disturbance  of  interest  so  that  the 
exact  form  of  the  conditions  will  not  be  critical  to  the  performance  of 
the  model.  For  ease  of  implementation  a  periodic  boundary  condition  is 
used.  That  is,  each  edge  becomes  a  neighbor  to  the  other  for  the 
computation  of  derivatives  and  mass  flux  terms. 


To  implement  the  periodic  boundary  condition  and  to  provide 
sufficient  grid  overlap  for  the  determination  of  the  depth  derivatives 
(see  Appendix  II),  the  bottom  three  grid  lines  are  repeated  at  the  top 
of  the  grid.  In  the  notation  of  Fig.  3,  the  grid  points  along  the  grid 
line  at  N-2  are  the  same  as  those  along  the  line  at  1.  The  lines  at  2 
and  N-l,  and  at  3  and  N  also  form  equivalent  pairs.  The  front  points 
are  also  subject  to  the  periodic  boundary  condition.  This  is  important 
primarily  for  smoothing.  When  a  front  point  moves  above  a  limit  at  the 
'top'  edge,  it  is  moved  to  the  'bottom'  of  the  grid,  as  though  it  were 
just  entering  the  grid  from  the  bottom. 

2.9  -  Parent  Pool  Boundary  Condition 

The  parent  pool  region  is  assumed  to  be  undisturbed  by  the 
motions  of  the  front.  The  depth  is  constant  at  DQ,  the  scaling  depth, 
so  that  5*1.  The  v  velocity  is  zero  since  the  parent  pool  and  ambient 
waters  are  assumed  to  have  no  relative  motion.  The  u  velocity  component 
compensates  (approximately)  for  the  mass  entrainment  in  the  frontal 
region.  Derivatives  of  all  of  the  variables  are  assumed  to  vanish  in 
the  parent  pool. 

These  conditions  are  applied  to  the  model  boundary  along  the 
parent  pool  as  initial  conditions.  However,  they  cannot  be  used  as 
fixed  boundary  conditions  at  later  time  steps.  If  the  boundary  values 
are  held  at  fixed  levels,  disturbances  which  are  generated  in  the  interior 
of  the  grid  will  not  be  able  to  escape  the  domain,  but  will  be  reflected 
or  trapped  by  the  boundary.  A  boundary  condition  which  allows  distur- 
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bances  to  propagate  freely  out  of  the  domain  is  required.  As  discussed 
in  section  2.8,  this  problem  is  solved  for  the  end  boundaries  through 
the  use  of  a  periodic  boundary  condition.  For  the  parent  pool  boundary, 
however,  there  is  no  other  boundary  available  with  which  to  form  the 
periodic  condition.  A  Sommerfeld  radiation  condition  is  used  instead. 
This  condition  may  be  expressed: 

“t  +  cux  *  0 

where  u  is  any  variable  and  c  is  the  phase  velocity  of  the  waves. 

Orlanski  (1976)  suggests  an  algorithm  for  this  condition  which 
is  employed  here.  This  method  consists  of  two  steps.  First,  for  each 
variable  a  disturbance  phase  speed  is  calculated  at  each  boundary  point 
from  the  neighboring  points  in  the  grid  interior.  The  values  of  the 
variables  along  the  boundary  are  then  extrapolated  from  the  neighboring 
points,  using  the  calculated  phase  speeds.  If  the  calculated  phase 
speed  exceeds  the  maximum  internal  wave  speed,  it  is  reset  to  equal  the 
internal  wave  speed.  Similarly,  negative  calculated  phase  speeds  (dis¬ 
turbances  moving  away  from  the  boundary)  are  reset  to  zero. 


THE  LAGRANGIAN  MODEL 


3.1  -  Basic  Equations 

Since  the  details  of  the  flow  in  the  turbulent  inner  region  are 
not  known,  the  movement  of  the  surface  expression  of  the  front  cannot  be 
determined  directly.  We  may  infer  the  frontal  position  from  the  position 
of  the  boundary  between  the  inner  and  outer  regions  by  assuming  that  the 
configuration  of  the  inner  region  in  the  cross  frontal  direction  is  not 
affected  by  the  movement  of  the  boundary.  Since  one  of  the  major  ob¬ 
jectives  of  this  research  is  to  follow  the  time  evolution  of  the  front, 
modeling  the  inner/outer  region  boundary  is  of  great  importance.  Dif¬ 
ferential  equations  may  be  written  describing  the  movement  of  the  free 
boundary  and  a  differencing  scheme  must  be  applied  in  order  to  integrate 
the  equations  through  time.  Since  an  Eulerian  approach  does  not  work 
well  for  free  boundaries,  we  link  a  Lagrangian  model  for  this  boundary 
to  the  Eulerian  model  for  the  outer  region. 

The  position  of  the  inner/outer  boundary  is  marked  by  a  series 
of  quasi-Lagrangian  drift  points.  These  points  are  initially  set 
slightly  closer  together  than  the  spacing  of  the  fixed  grid.  The 
initial  conditions  of  frontal  disturbance  wavelength  and  amplitude  also 
determine  the  initial  positions  of  the  points  along  the  grid.  These 
drift  points  are  Lagrangian  in  that  their  velocities  are  determined  from 
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the  Lagrangian  flow  equations.  However,  their  motion  is  determined  not 
only  by  the  Lagrangian  velocities,  but  by  the  mass  flux  across  the 
boundary  as  shown  in  Fig.  4.  The  drift  points  move  with  the  boundary, 
defining  its  position.  They  do  not  stay  with  the  initial  water  parcels 
as  fully  Lagrangian  drifters  would,  except  in  the  case  of  zero  mass  flux 
across  the  boundary. 

Along  the  inner/outer  boundary  the  depth  of  the  interface  is 
held  constant,  0  =  0^.  Non-dimensionally  this  is 


A  fluid  particle  which  is  initially  on  the  boundary  will  move  as 
a  Lagrangian  drifter  in  the  velocity  field.  That  is,  the  displacement 
through  time  will  be 

=  up(xp,JV^  (3.1.2a) 

Ht^p^))  *  (3.1.2b) 

3.2  -  Mass  Flux  Correction 

If  the  fluid  parcel  were  to  remain  on  the  line  defining  the 
boundary,  the  equations  (3.1.2)  would  describe  the  displacement  of  the 
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boundary.  However,  turbulent  entrainment  may  occur  in  the  inner  region, 
resulting  in  a  net  gain  (or  loss)  of  water  in  that  region.  If  the 
volume  of  water  in  the  inner  region  is  to  remain  constant,  then  the 
inner/outer  boundary  must  move  with  respect  to  the  local  water  at  a 
velocity  which  compensates  for  the  entrainment.  If,  in  a  unit  of  time  a 
volume  V  is  entrained  into  a  section  of  the  inner  region  of  unit  width, 
then  the  velocity  of  the  boundary  relative  to  the  water  must  be 


V 


r 


(3.2.1) 


We  note  that  there  may  also  be  a  velocity  shear  applied  to  the 

boundary.  This  shear  velocity,  v  ,  will  displace  water  parcels  from  a 

point  fixed  on  the  boundary  in  the  same  manner  as  the  entrainment 
velocity,  \>r.  Since  the  boundary  may  be  curved  and  not  aligned  with  the 

model  grid,  a  transformation  is  required  to  relate  vr  and  to  the 

A  A 

velocities,  u^  and  v^,  which  are  aligned  with  the  grid  and  are  determined 


uf  =  vr  cos  9n  -  s1n  sn 

(3.2.2a) 

vf  =  vr  sin  *  vj  cos  9„ 

(3.2.2b) 

This  relation  is  shown  graphically  in  Fig.  5.  The  instantaneous 
velocity  of  a  point  which  is  fixed  on  the  boundary  will  be  equal  to  the 
velocity  of  a  water  parcel  on  the  boundary  plus  the  velocity  of  the 
boundary  with  respect  to  the  water. 
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ub  =  u  ”  uf  (3.2.3a) 

vb  =  v  “  vf  (3.2.3b) 

If  small  displacements  are  assumed,  this  Eulerian  scheme  may  be 
linked  to  a  Lagrangian  scheme  to  give  the  displacements  in  time  of  the 
boundary  points 


ft(xb(t))  =  ‘  Uf  (3.2.4a) 

=  vp^xp»'yp,t)  "  vf  (3.2.4b) 


In  the  small  time  increments  involved  in  the  finite  difference 
method  of  the  model,  the  small  displacement  condition  can  easily  be  met. 
Since  the  depth  field  is  available  as  a  boundary  condition  at  the 
boundary  and  as  a  model  result  on  the  model  grid  in  the  outer  region  and 
from  extrapolations  into  the  inner  region  from  the  boundary  (see  Appendix 
I),  the  Up  and  vp  values  may  be  obtained  from  the  momentum  equations 


u 


v 


pr 

+  U  U  +  V  u 

P  Pn  P  PC 

+  U  V 

> 

> 

4- 

PT 

P  Pn 

P  PC 

+  u 

P 


(3.2.3a) 

(3.2.5b) 


or, 


dt 


(3.2.6a) 


u  +  u 
0  00 


(3.2.6b) 


3.3  -  Difference  Equations 


To  fit  the  general  differencing  cycle  of  the  overall  model,  the 

momentum  and  displacement  equations  must  be  put  into  a  difference  form. 

The  difference  forms  of  (3.2.6)  are 

Up(T  +  At)  =  Up(t)  +  Ar(Vp  -  V,  - 

(3.3.1a) 

Vp(T  +  At)  *  Vp ( t )  +  AT(-Up  +  U^  -  5?) 

(3.3.1b) 

This  gives  the  Lagrangian  velocity  of  the  water  along 

boundary.  The  motion  of  the  drift  points  from  (3.2.4)  is 

the 

Xjj(t  +  At)  =  xb(T)  +  At ( up  -  uf) 

(3.3.2a) 

yb(t  +  At)  =  yb(t)  +  Ar(vp  -  vf) 

(3.3.2b) 

This  formula,  applied  to  each  point  gives  its  position  on  the 
fixed  grid  at  each  time  step.  These  difference  equations  are  solved  by 
a  predictor-corrector  scheme.  The  time  step  used  is  the  same  as  that 
used  in  the  finite  differencing  scheme  for  the  fixed  grid,  so  that  the 
two  solutions  can  advance  stepwise  in  parallel  through  time. 


INITIAL  CONDITIONS 


4.1  -  Depth 

Following  Formato  (1979),  we  may  derive  the  form  of  the  depth 
profile  from  (2.3.3),  (2.7.3),  and  (2.3.6).  To  find  the  basic  equilibrium 
state  we  first  eliminate  time  derivatives.  Then,  taking  y  derivatives 
to  be  small  when  compared  to  x  derivatives,  we  set  u,  =  v  *  0.  The  x 

s  t 

momentum,  vorticity,  and  continuity  equations  are  then 


uu 

! 

-  v  + 

6  3  0 
n 

(4.1.1) 

fv  +  1' 

_ 0 _ 

5 

3  0 

(4.1.2) 

1 

n 

(U6)n 

*  0 

(4.1.3) 

Equations  (4.1.2)  and 

(4.1.3; 

)  may  be 

integrated  to  yield 

V  +  1 

n 

3  6 

(4.1.4) 

u6  3 

uo 

(4.1.5) 

Differentiating  (4.1.1)  with  respect  to  n  and  substituting  for 
v^  from  (4.1.4)  and  for  u  from  (4.1.5)  gives  an  equation  for  5  alone. 
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(4.1.6) 


An  exact  solution  may  be  obtained  for  this  equation,  but  the 

solution  is  transcendental .  However,  by  noting  that  uQ  =  F  ,  the  cross 

2 

stream  Froude  number  and  requiring  Fp  <<  5Q  <  1 ,  we  may  form  a  small 
perturbation  expansion  of  the  form 

6(n)  *  5^(n)  +  F25^(n)  +  •••  (4.1.7) 

Applying  (4.1.7)  to  (4.1.6),  we  have  for  the  lowest  order 

solution 


5(0)  -  +1*0  (4.1.3) 

nn 

We  introduce  a  new  coordinate  in  the  x  direction  n$  =  (n  - 
nb)/RQ,  where  nb  is  the  location  of  the  inner/outer  boundary  and  Rq  is 
the  baroclinic  Rossby  radius.  At  the  boundary  ns  -  0.  The  boundary 
conditions  to  be  met  are 


S(0)(ns  *  0)  -  Sb  (4.1.9a) 

lim  5(0)(n  )  =  1  (4.1.9b) 

n$  ♦  » 

Subject  to  these  conditions  the  solution  of  (4.1.8)  is 
o(n$)  *  1  -  (1  -  5b)e  s 


(4.1.10) 
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This  gives  depth  as  a  function  of  x  directed  distance  from  the 
front.  It  is  the  exponential  relation  derived  by  Stommel  (1976). 

4.2  -  Velocity 

An  expression  for  the  u  component  of  velocity  may  be  obtained 
directly  from  (4.1.5) 


This  expresses  cross  stream  conservation  of  mass.  At  the  parent 
pool  boundary  u  =  uQ  =  Fr,  and  at  the  free  boundary  u  =  t»0/6b . 

2 

In  (4.1.1)  the  first  term  will  be  small  by  a  factor  of  Fr  so 
that  it  may  be  neglected  to  lowest  order.  The  v  component  of  velocity 
from  this  is 

v  =  5  (4.2.2) 

n 

The  initial  conditions  thus  include  geostrophic  balance  for  the 
v  component.  Using  the  relation  for  the  depth  given  above,  (4.2.2)  may 
be  expressed  in  a  slightly  different  form.  Differentiating  (4.1.10) 

«  -  (1  -  3b)e  ^  (4.2.3) 

or 

<5  =l-o 
n 


(4.2.4) 
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so  that 

v  =  1  -  <5  (4.2.5) 

This  relation  will  be  useful  later  in  deriving  a  smoothing 
scheme  (Appendix  I). 


MODEL  TESTING 


5.1  -  General  Description 

The  model  described  above  and  in  the  Appendices  was  coded  in 
FORTRAN  and  run  on  a  DEC-10  computer  at  the  University  of  Delaware. 

After  a  set  of  parameters  was  chosen,  the  initial  conditions  were 
computed  and  the  model  was  run  for  one  time  step  with  the  free  boundary 
held  in  a  fixed  position  and  the  solutions  in  the  inner  region  held  at 
the  initial  values.  This  was  done  to  allow  for  any  minor  adjustments  to 
the  initial  conditions  necessary  in  the  outer  region.  Following  this 
initialization,  the  model  runs  were  executed  piecewise;  that  is,  the 
model  was  run  for  one  or  more  inertial  periods,  the  results  were  exam¬ 
ined,  and  the  execution  was  resumed,  continuing  for  one  (or  more) 
inertial  periods(s)  from  the  point  at  which  it  had  been  halted.  This 
pattern  was  repeated  until  a  sufficient  total  time  period  had  been 
modeled. 


The  front  was  established  parallel  to  the  y  axis,  with  the 
perturbations  in  the  x  direction  as  shown  in  Fig.  2.  The  model  grid  was 
set  up  with  one  hundred  points  in  the  y  direction  (with  three  points 
overlapping  as  described  in  section  2.8),  and  forty  points  in  the  x 
direction.  To  specify  a  location  we  may  use  a  grid  coordinate  system 
based  on  these  points.  In  this  system  the  x  coordinate  ranges  from 
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one  to  forty  and  the  y  coordinate  ranges  from  one  to  one  hundred,  with 
the  point  (1,1)  located  in  the  lower  left  corner.  For  most  of  the  cases 
run,  the  unperturbed  portion  of  the  front  was  initially  set  at  an  x 
position  of  3.5  in  grid  coordinates.  For  some  cases  the  initial  frontal 
position  was  set  at  a  larger  x  value  to  allow  room  in  the  grid  for  the 
front  to  move  to  the  left.  This  condition  arose  when  the  frontal 
propagation  speed  was  negative. 

The  testing  of  the  model  was  performed  with  a  set  of  parameter 
values  characteristic  of  the  Gulf  Stream  downstream  of  its  point  of 
separation  from  the  continental  margin.  The  parameter  values  used  for 
the  testing  are  shown  in  Table  1. 

This  choice  of  parameters  determined  values  for  several 
related  variables.  These  are  shown  in  Table  2. 

The  disturbance  initially  imposed  on  the  front  had  a  wave¬ 
length  of  300  km  and  an  amplitude  of  100  km.  The  cross-stream  Froude 
number  was  set  at  -0.02.  This  value  was  chosen  because  preliminary 
tests  indicated  that  the  model  was  most  unstable  near  the  inertial 
frequency  with  a  negative  Froude  number,  and  -0.02  was  the  largest 
value  that  seemed  physically  reasonable. 

Several  types  of  stability  must  be  considered  here.  The  first 
is  the  numerical  stability  of  the  numerical  scheme.  An  inaccurate 
numerical  scheme  can  yield  incorrect  results  in  the  form  of  a  divergent 
solution  ('blow  up')  which  is  easily  detectable.  No  errors  of  this  type 
occurred  for  any  run.  A  numerical  scheme  can  also  give  an  erroneous 


38 


Table  1 

Model  parameters  for  the  Gulf  Stream  front 


Grid  Spacing  (x  and  y) 
Time  Step 

Stratification  (Ap/p^) 
Latitude 


10  km 
855.974  s 
0.0013 
39°  N 
300  m 


900  m 
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Table  2 

Physical  characteristics  of  the  modeled  frontal  region 


Inertial  period  19.02  hr 

2 

g'  (reduced  gravity)  0.0127  m/s 

c  (maximum  internal  wave  speed)  3.386  m/s 

Maximum  internal  wave  speed  at  the  front  1.955  m/s 

x  (baroclinic  Rossby  radius)  36.9  km 
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result  in  the  fora  of  a  non-convergent  solution  which  is  more  difficult 
to  detect.  The  model  was  checked  for  non-convergence  by  testing  the 
dependence  of  the  results  on  the  size  of  the  time  step  and  the  grid 
spacing,  as  discussed  in  section  5.3. 

The  second  stability  consideration  is  the  response  of  the 
model  to  oscillations  near  the  inertial  frequency.  Numerical  models 
coiranonly  respond  strongly  to  excitation  near  'natural'  frequencies, 
such  as  the  inertial  frequency.  For  example,  in  a  general  circulation 
model  investigating  baroclinic  instability,  Orlanski  and  Cox  (1973) 
observed  significant  inertial  frequency  oscillations  along  the  front 
associated  with  the  Florida  Current.  If  there  is  no  mechanism  in  the 
model  to  dissipate  the  energy  in  these  oscillations,  as  there  must  be 
in  the  ocean,  they  can  grow  with  time  and  eventually  dominate  all  cf 
the  motion  in  the  model.  Since  the  main  interest  of  this  investigation 
involves  oscillations  at  frequencies  much  lower  than  inertial,  oscil¬ 
lations  near  the  inertial  frequency  are  of  secondary  importance.  To 
keep  the  high  frequency  motions  from  subordinating  the  low  frequency 
motions,  the  high  frequency  motions  must  be  suppressed.  This  is  done 
through  the  use  of  interfacial  friction,  as  discussed  in  section  5.2. 
Ideally,  we  would  wish  to  restrict  the  inertial  oscillations  from 
unlimited  growth  without  suppressing  them  entirely. 

The  third  type  of  stability  directly  pertains  to  the  frontal 
perturbations  which  are  the  subject  of  this  research.  These  pertur¬ 
bations  are  spatially  large  (tens  of  km)  and  temporally  slow  (several 
inertial  periods).  The  main  question  to  be  asked  here  is:  under  what 
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circumstances  are  these  perturbations  unstable?  The  criteria  of  sta¬ 
bility  may  be  set  in  two  ways.  Since  the  initial  amplitude  of  an 
applied  perturbation  is  known,  the  amplitude  at  a  later  time  may  be 
compared  to  see  if  it  is  increasing  (unstable)  or  decreasing  (stable). 
Alternately,  since  the  initial  kinetic  and  potential  energies  are  known, 
the  ratios  of  perturbation  energy  to  mean  energy  can  be  calculated  at  a 
later  time  to  see  if  they  are  increasing  (unstable)  or  decreasing 
(stable) . 


Imposing  the  latter  two  criteria  for  stability  can  (and  some¬ 
times  will)  lead  to  conflicting  conclusions.  We  may  relate  this  con¬ 
flict  by  analogy  to  a  surface  gravity  wave  approaching  a  shore  over  a 
bottom  with  a  shallow  slope.  In  that  system  kinetic  energy  is  lost  due 
to  turbulence  within  the  fluid  and  friction  along  the  bottom.  Although 
these  losses  are  slow,  they  occur  most  rapidly  where  the  deviations  from 
the  mean  are  the  greatest,  i.e.  near  the  wave  crest.  Meanwhile,  the 
amplitude  of  the  wave  increases  with  time.  In  the  case  of  a  surface 
gravity  wave  all  of  the  energy  of  the  system  is  wave  (fluctuation) 
energy,  so  that  fluctuation  to  mean  energy  ratios  are  not  meaningful. 
However,  it  can  be  seen  that  the  total  fluctuation  energy  will  decrease 
over  time,  indicating  stability.  At  the  same  time,  however,  the  surface 
slope  and  the  amplitude  of  the  wave  are  increasing,  indicating  instabil¬ 
ity.  Thus,  the  stability  of  the  interface,  exhibited  by  the  amplitude, 
does  not  necessarily  follow  the  stability  of  the  interior  flow,  ex¬ 
hibited  by  the  kinetic  energy. 
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In  this  study  the  displacement  amplitude  criterion  is  used 
because  it  relates  to  other  properties  that  are  computed,  e.g.,  phase 
speed,  and  also  relates  to  observable  properties.  The  energy  criterion 
is  used  to  examine  the  dynamics  of  the  flow  and  to  allow  broader  con¬ 
clusions  to  be  drawn  regarding  the  nature  of  the  instabilities. 

5.2  -  Selection  of  the  Coefficient  of  Friction 

As  described  above,  interfacial  friction  is  used  in  this  model, 
in  part,  to  stabilize  the  computations  with  respect  to  inertial  oscil¬ 
lations.  The  corresponding  coefficient  of  friction,  Cq,  cannot  be 
determined  readily  from  observational  results.  Several  values  were 
tried  in  the  model  to  determine  the  value  best  suited  for  damping  the 
spurious  inertial  oscillations.  Tests  were  run  for  friction  coeffi¬ 
cients  of  0.0010,  0.0015,  0.0020,  and  0.0040.  Plots  of  the  corresponding 
kinetic  energy  fluctuations  with  time  are  shown  in  Fig.  oa-d.  These 
plots  express  the  kinetic  energy  as  normalized  by  the  initial  kinetic 
energy  lv/£v  .  The  variation  in  kinetic  energy  shows  that  the  two 
smaller  values  of  Cg  led  eventually  to  increasing  inertial  frequency 
oscillations  which  led  rapidly  to  velocities  in  excess  of  c  and  thus  to 
termination  of  the  computations,  while  the  largest  value  showed  over¬ 
damping.  Accordingly,  a  value  of  0.002  was  selected  for  Crj  for  all  of 
the  model  runs. 

The  use  of  interfacial  friction  makes  the  bottom  boundary  of  the 
model  (the  interface)  a  net  sink  of  kinetic  energy.  Although  some 
energy  can  be  introduced  into  the  model  at  the  parent  pool  boundary  in 


43 


the  form  of  potential  energy,  this  was  not  enough  in  these  tests  to 
compensate  for  the  energy  loss.  The  percentage  losses  of  kinetic  and 
potential  energy  as  a  function  of  Cq  are  shown  in  Fig.  7.  The  x 
axis  intercepts  indicate  that  at  a  small  CQ  the  net  losses  would 
vanish.  However,  a  value  of  Cq  that  small  would  not  have  been  effec¬ 
tive  in  suppressing  the  inertial  oscillations. 

Net  losses  of  both  kinetic  and  potential  energy  were  seen 
for  all  cases.  The  differences  in  the  rates  of  loss  between  cases 
was  not  caused  primarily  by  differences  in  the  evolution  of  the 
frontal  perturbations.  Therefore,  the  absolute  magnitudes  of  the 
energies  were  not  used  for  comparison.  Instead,  ratios  of  the 
fluctuation  to  mean  energies  (as  determined  by  Eq.  2. 6. 3-2. 6.6)  were 
used  to  compare  the  staoility  characteristics.  It  should  be  noted 
that  the  fluctuations  measured  thus  are  with  respect  to  means  taken 
in  the  y  direction, 

5.3  -  Sensitivity  Tests 

As  stated  above,  the  standard  time  step  and  grid  spacing  were 
set  at  855.974  seconds  and  10.  km,  respectively.  To  ensure  that  the 

* 

results  of  the  model  were  convergent  and  not  sensitive  to  the  choice 
of  either  of  these,  test  runs  were  made  with  different  values  for 
these  parameters.  The  maximum  time  step  allowed  for  numerical  stability 
is  fixed  by  the  grid  size  and  the  speed  of  the  fastest  internal  wave. 

The  ratio  of  the  time  step  used  to  this  maximum  time  step  is  the 
Courant  number,  and  must  be  less  than  or  equal  to  one.  For  the 
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standard  case  the  Courant  number  is  0.82.  To  test  the  time  step 
sensitivity  a  run  was  made  with  a  Courant  number  of  0.65,  or  a  time 
step  of  570.649  seconds.  For  the  grid  spacing  sensitivity  test  an 
8  km  grid  was  used.  The  Courant  number  for  this  latter  test  was  set 
at  0.82,  as  in  the  standard  case,  so  the  time  step  was  then  684.9 
seconds. 


The  three  cases,  one  standard  and  two  comparison  tests,  were 
run  for  ten  inertial  periods.  This  required  800  time  steps  for  the 
standard  case,  1000  for  the  grid  sensitivity  test,  and  1200  for  the 
time  step  sensitivity  test.  Contours  of  depth  at  the  initial  time 
and,  for  each  test  case,  after  ten  inertial  periods  are  shown  in  plan 
view  in  Fig.  8.  The  major  features  of  these  three  figures  are  quite 
similar,  although  there  are  small  differences.  The  two  cases  with  the 
higher  Courant  number  have  a  greater  amount  of  short  wavelength 
variation  in  the  deep  region  as  shown  by  the  oscillations  of  the  0.9 
contour  line.  Along  the  leading  side  of  the  disturbance  the  low 
Courant  number  case  appears  to  be  slightly  more  rounded.  Minor 
variations  are  also  seen  elsewhere  along  the  front. 

Other  physical  properties  of  the  results  were  computed  for 
intercomparison.  The  phase  speed  was  determined  by  following  the 
movement  of  the  point  of  maximum  disturbance.  The  position  of  this 
point  was  determined  from  a  least  squares  fit  of  a  cubic  of  the  five 
points  with  the  largest  x  coordinate  values.  The  phase  speed  was 
computed  from  a  linear  least  squares  fit  of  a  series  of  these  points 
sampled  every  half  inertial  period.  Positive  values  mean  that  the 


45 


wave  is  moving  towards  the  top  (downstream)  of  the  grid.  The  frontal 
propagation  speed  was  taken  as  the  speed  in  the  x  direction  of  the 
straight,  unperturbed  portion  of  the  front  near  the  top  and  bottom 
boundaries.  Positive  values  mean  that  the  front  is  moving  towards 
the  right  (parent  pool).  The  amplitude  was  calculated  as  the  differ¬ 
ence  in  x  position  between  the  point  of  maximum  disturbance  (crest) 
and  a  point  determined  similarly  as  the  smallest  x  position  point 
along  the  front  (trough).  These  results  are  compared  in  Table  3. 

These  figures  are  quite  close,  generally  within  2%,  and 
indicate  that  the  dependence  of  the  model  on  the  differencing  parameters 
is  slight.  The  grid  sensitivity  test  showed  the  greatest  variation, 
indicating  that  greater  spatial  resolution  would  be  a  more  important 
improvement  than  a  decreased  Courant  number.  Greater  resolution  would 
also  suppress  the  short  wavelength  variations  shown  in  Figs.  6b)  and 
c). 


The  ratios  of  fluctuation  energy  to  mean  energy  were  computed 
for  both  the  kinetic  and  potential  energies.  The  initial  and  final 
(after  ten  inertial  periods)  values  of  these  ratios  for  the  three  test 
cases  are  shown  in  Table  4. 

The  results  for  the  standard  and  the  time  sensitivity  cases 
were  again  quite  close,  within  4%.  The  values  of  the  ratios  for  the 
grid  sensitivity  test  were  uniformly  higher  than  the  corresponding 
values  for  the  other  two  cases  since  the  grid  sensitivity  test  used 
the  same  number  of  grid  cells  as  the  other  tests  but  a  smaller  grid 
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Table  3 

Disturbance  propagation  results  for  sensitivity  tests 


Standard 

Time 

Sensitivity 

Grid 

Sensitivity 

Phase  speed  (cm/ sec) 

3.94 

3.96 

3.89 

Frontal  speed  (cm/sec) 

2.42 

2.27 

2.37 

Amplitude  (km) 

108.3 

106.8 

106.5 
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Table  4 


Fluctuation  to 

mean  energy 

ratios  for  sensitivi 

ty  tests 

Standard 

Time 

Sensitivity 

Grid 

Sensitivity 

Kinetic  (xlO-^) 

Initial 

119. 

119. 

147. 

Final 

336. 

351. 

443. 

%  Increase 

282. 

295. 

301. 

Potential  (xlCf4) 

Initial 

29.2 

29.2 

32.3 

Final 

53.1 

51.7 

65.9 

"  Increase 

182. 

177. 

204. 
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spacing,  the  total  horizontal  span  of  the  model  was  reduced.  The 
reduced  span  caused  the  applied  perturbation  to  occupy  a  larger  fraction 
of  the  model  grid.  Therefore,  the  perturbation  components  of  the 
energies  were  a  larger  fraction  of  the  total  energies.  Proportionally, 
the  differences  between  the  initial  and  final  ratios  for  the  grid 
sensitivity  case  were  quite  similar  to  the  ratios  for  the  other  two 
cases . 


These  sensitivity  tests  demonstrate  that  the  model  is 
numerically  convergent.  The  set  of  differencing  parameters  used  is 
acceptable  and  the  results  were  not  seriously  affected  by  the  choice 
of  either  the  Courant  number  or  the  grid  spacing.  Thus,  for  the  model 
runs  intended  to  approximate  the  conditions  of  the  Gulf  Stream  the 
Table  1  values  were  used. 


MODEL  RESULTS 


6.1  -  Introduction 

With  the  model  verified  as  numerically  convergent  and  with 
the  growth  of  the  high  frequency  motions  controlled  by  friction,  the 
effect  of  variations  in  the  free  parameters  on  the  stability  of  the 
flow  may  be  examined.  In  these  tests  conditions  characteristic  of  two 
regions  where  persistent  fronts  occur  were  used.  The  regions  were  the 
Gulf  Stream  downstream  of  its  point  of  separation  from  the  continental 
margin,  and  the  Sargasso  Sea  at  30°  N.  These  regions  were  defined  by 
the  latitude,  parent  pool  depth,  inner/outer  boundary  depth  and 
stratification.  These  regions  were  chosen  to  represent  two  different 
oceanographic  regimes.  The  Gulf  Stream  is  a  swift,  deep  current  which 
is  the  predominant  dynamic  feature  in  its  region.  The  currents  con¬ 
nected  with  the  Sargasso  Sea  fronts  are  much  slower  and  shallower. 

They  appear  to  have  a  seasonal  existence  as  opposed  to  the  permanence 
of  the  Gulf  Stream  (Voorhis,  1969).  They  are,  however,  sufficiently 
long  lived  (several  months)  to  qualify  as  permanent  for  the  purposes 
of  this  discussion. 

In  modelling  fronts  in  these  two  regions  several  other 
parameters  were  expected  to  influence  the  flow.  For  both  regions  the 
cross  stream  Froude  number  ana  the  amplitude  of  the  initial  disturbance 
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were  varied  to  determine  the  effect  on  the  propagation  of  the  distur¬ 
bance.  For  the  Sargasso  Sea  front  the  wavelength  of  the  initial 
disturbance  was  also  varied. 

6.2  -  Gulf  Stream  Unperturbed  Front 

Before  testing  the  stability  characteristics  of  a  perturbed 
front,  the  general  characteristics  of  an  unperturbed  (straight)  front 
were  examined.  The  characteristics  of  greatest  interest  were  the 
frontal  cross  stream  propagation  speed  and  the  change  in  the  kinetic 
and  potential  energies  with  time.  Trials  were  made  for  three  cases 
with  varying  cross  stream  Froude  numbers,  Fr  *  -0.02,  Fp  =  0.00,  Fr  =  0.01. 
The  F  *  0.00  value  was  chosen  to  provide  a  test  with  no  induced  cross 
flow  so  that  the  'basic'  dynamics  could  be  studied.  A  Froude  number 
with  a  magnitude  of  0.02  was  considered  the  largest  that  was  physically 
realistic.  Both  positive  and  negative  Froude  numbers  were  used  to  test 
the  effect  of  the  direction  of  the  cross  stream  flow,  but  a  value  of 
0.02  yielded  a  frontal  propagation  speed  that  was  too  large  to  be 
plausible  (see  below).  Accordingly,  the  positive  Froude  number  was 
reduced  to  0.01,  while  the  negative  remained  at  -0.02. 

For  each  of  the  three  Froude  number  cases  the  model  was  run 
for  ten  inertial  periods.  The  trials  with  Fr  =  -0.02  and  Fr  *  0.00  were 
started  with  the  front  at  ar.  x  position  of  3.5  grid  coordinates.  The 
trial  with  Ff  =  0.01  was  run  with  the  front  initially  at  an  x  position 
of  5.5  grid  coordinates.  To  test  the  (possible)  effects  of  the  model 
boundaries  on  the  frontal  movement  an  extra  run  was  made  for  Fr  =  0.01, 
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started  with  the  front  at  an  x  position  of  3.5  grid  coordinates.  This 
test  was  run  for  four  inertial  periods.  The  left  boundary  cannot  affect 
the  results  of  this  model  as  long  as  the  front  is  at  least  one  grid 
cell  away  from  it.  The  right  (parent  pool)  boundary  might  affect  the 
flow  if  it  is  not  a  sufficiently  large  distance  from  the  front.  A 
significant  difference  between  the  two  *  0.01  cases  would  indicate 
an  insufficiently  large  model  domain.  The  position  of  the  front  as  a 
function  of  time  for  each  of  these  four  cases  is  plotted  with  open 
symbols  in  Fig.  9.  The  four  lines  with  solid  symbols  at  the  end  points 
are  least  squares  linear  fits  to  the  frontal  positions,  taken  once  per 
inertial  period. 

For  the  negative  Frouae  number  case  the  front  moved  toward 
the  right  (larger  x  values).  For  the  other  two  cases  the  motion  was 
to  the  left  (smaller  x  values).  For  all  of  the  cases  the  motion  of 
the  front  was  nearly  uniform  throughout  the  run.  Frontal  propagation 
speeds  were  calculated  from  the  least  squares  fits  and  are  shown  in 
Table  5.  This  table  also  includes  the  results  of  a  test  with  Fr  =  -0.02 
for  which  the  initial  position  of  the  front  was  along  x  =  12.5  grid 
coordinates. 

The  difference  in  initial  position  for  the  Fp  *  0.01  cases  had 
a  negligible  effect  on  the  frontal  propagation  speeds.  The  greater 
difference  in  starting  position  between  the  two  Fr  s  -0.02  cases  affected 
the  speed  less  than  3%.  This  result  indicates  that  the  model  domain  is 
large  enough  in  the  x  direction  so  that  the  presence  of  the  parent  pool 
boundary  does  not  significantly  affect  the  motion  of  the  free  boundary. 
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Table  5 

Frontal  propagation  speeds  for  an  unperturbed  front 


Froude  number 

Initial  position 

Speed  (cm/s) 

-0.02 

3.5 

2.41 

-0.02 

12.5 

2.34 

0.00 

3.5 

-2.48 

0.01 

3.5 

-4.94 

0.01 


5.5 


-4.98 
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Thus,  the  initial  x  position  was  set  at  3.5  for  the  tests  in  which  the 
front  moved  to  the  right.  For  those  tests  in  which  the  front  moved  to 
the  left,  the  initial  x  position  was  set  at  a  sufficiently  large  value 
to  keep  the  front  in  the  grid  for  the  entire  test. 

The  zero  Froude  number  case  indicates  that  the  flow  was  not 
exactly  in  geostrophic  balance.  The  deviation  from  geostrophy  was 
probably  caused  in  part  by  the  interfacial  friction  which  slightly 
reduced  the  downstream  velocity.  If  the  downstream  velocity  is  not 
sufficient  to  geostrophically  balance  the  pressure  gradient  from  the 
slope  of  the  interfacial  depth,  a  compensating  cross  stream  flow  will 
develop.  Such  a  cross  stream  flow  was  observed  and  was  influenced  by 
the  Froude  number,  as  shown  in  Fig.  10.  The  velocities  shown  in  Fig.  10 
are  relative  to  the  front,  which  is  moving,  rather  than  to  the  fixed 
grid. 

The  results  in  Table  5  show  that  the  presence  of  basic  state 
across  frontal  flow  (non-zero  Froude  number)  affects  the  propagation  of 
the  front  across  the  grid.  The  changes  in  the  cross  stream  frontal 
propagation  speed  are  only  a  few  percent  of  the  downstream  velocity, 
which  indicates  that  they  are  second  order.  The  mechanisms  causing 
these  changes  from  the  zero  Froude  number  case  are  not  clear,  but  were 
probably  connected  with  the  variations  in  the  relative  cross  stream 
velocity  shown  In  Fig.  10.  Although  all  of  the  cases  showed  positive 
relative  cross  stream  velocities  and  volume  losses,  the  negative  Froude 
number  case  (dashed  line)  had  smaller  velocities  than  either  the 
Fr  3  0.00  or  Fr  =  0.01  (dotted  line)  cases.  In  spite  of  a  loss  in 


volume  and  a  positive  relative  cross  stream  velocity,  the  Fr  =  -0.02 
case  showed  a  rightward  movement  in  the  front. 

The  series  of  frontal  positions  shown  in  Fig.  9  shows  that  the 
flow  supporting  the  cross  frontal  propagation  of  the  front  developed 
rapidly  and  continued  uniformly  through  time.  In  the  range  of  F^ 
examined  this  effect  was  linear;  that  is,  the  change  in  the  frontal 
propagation  speed  was  proportional  to  the  cross  stream  Froude  number, 
approximately  2.5  cm/s  for  each  0.01  increment  in  the  Froude  number. 
Since  the  frontal  propagations  for  the  zero  and  negative  Froude  number 
cases  were  in  opposite  directions,  a  negative  Froude  number  may  be 
necessary  for  the  front  to  maintain  a  fixed  position  on  the  grid  over 
a  long  period  of  time.  These  results  suggest  a  value  of  -0.01  for  the 
Froude  number.  Kao  (1980)  obtained  an  equilibrium  state  for  the  Gulf 
Stream  with  a  cross  flow  of  approximately  the  same  magnitude.  The 
cross  stream  flow  was  required  for  the  maintenance  of  the  frontal 
volume.  He  suggested  that  an  alongfront  pressure  gradient  (which  could 
not  be  used  here  due  to  the  periodic  boundary  conditions)  would  serve 
to  hold  the  front  stationary. 

Another  important  result  was  the  change  in  the  kinetic  and 
potential  energies  with  time.  These  were  calculated  from  the  model 
results  using  (2.6.3)-(2.5.4)  at  the  initial  and  final  (after  ten 
inertial  periods)  time  steps.  As  shown  in  Fig.  10  the  downstream 
velocity,  and  hence  the  kinetic  energy,  is  greatest  near  the  front  and 
diminishes  to  a  negligible  value  at  the  parent  pool.  The  depth,  hence 
the  volume  and  potential  energy,  is  smallest  near  the  front  and  rises 
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to  a  near  constant  value  towards  the  parent  pool.  All  three  cases 

showed  some  flattening  of  the  interface  and  thus  some  loss  of  volume 

in  the  frontal  region.  Kinetic  energy  losses  occurred  in  all  cases 

due  to  the  influence  of  interfacial  friction,  as  discussed  in  section 

5.2.  The  losses  for  these  tests  were  within  3%  of  each  other.  Small 

losses  also  occurred  in  the  potential  energies  consistent  with  the 

decreases  in  volume.  Slight  variation  between  the  tests  occurred. 

The  F„  =  -0.02  test  showed  the  least  volume  decrease  and  the  F„  =  0.01 
r  r 

test  showed  the  most. 

6.3  -  Gulf  Stream  Small  Amplitude  Perturbation 

To  study  the  behavior  of  perturbations  on  the  Gulf  Stream 
front,  tests  were  run  with  two  values  for  the  amplitude  to  wavelength 
ratio  for  the  initial  disturbance.  The  first  set  of  tests  had  a  small 
amplitude  to  wavelength  ratio  of  0.0333,  with  an  amplitude  of  10  km 
and  a  wavelength  of  300  km.  These  tests  were  run  with  the  same  set  of 
cross  stream  Froude  numbers  that  were  used  in  the  unperturbed  front 
cases.  They  were  run  for  only  six  inertial  periods  due  to  the  develop¬ 
ment  of  occluding  instabilities  along  the  front,  as  will  be  discussed 
later. 


The  perturbations  initially  applied  to  the  free  boundary  were 
centered  lengthwise  along  the  front.  Thus,  portions  of  the  front  at 
the  top  and  bottom  edges  of  the  model  grid  were  initially  straight  and 
undisturbed  by  the  frontal  perturbation.  This  was  done  to  eliminate 
edge  effects  near  the  perturbation  and  to  allow  room  for  movement  of 
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the  disturbance  along  the  front.  These  unperturbed  portions  of  the 
front  were  also  used  to  measure  the  frontal  propagation  speed  in  the 
cross  frontal  direction  for  comparison  with  the  unperturbed  front 
results.  Table  6  shows  the  initial  x  positions  of  the  front  in  grid 
coordinates  and  the  resulting  frontal  propagation  speeds. 

Different  initial  positions  were  used  to  allow  room  for  the 
leftward  propagating  fronts.  These  results  agree  very  well  with  those 
in  Table  5.  This  indicates  that  the  dynamics  of  the  regions  well  away 
from  the  perturbation  are  not  affected  by  a  perturbation  of  this  scale. 

Time  sequences  of  the  frontal  position  are  shown  in  Fig.  11  a-c) 
for  Fr  =  -0.02,  0.00,  and  0.01,  respectively.  Since  the  perturbation  is 
very  small,  the  x  coordinate  in  all  cases  has  been  expanded  by  a  factor 
of  five.  The  time  sequence  starts  with  the  leftmost  curve  of  each 
figure.  The  frontal  position  is  displayed  once  per  inertial  period 
thereafter.  These  curves  show  similar  patterns  of  deformations  at  the 
back  of  the  bulge  (toward  the  bottom  side  in  the  figure)  evolving  into 
short,  steep  disturbances  which  propagate  very  rapidly  to  the  leading 
side  of  the  bulge.  This  tendency  was  very  mild  for  Fr  =  -0.02  compared 
with  the  other  two  cases.  Small  disturbances  also  formed  downstream 
of  the  bulge.  In  two  cases  (F  *  -0.02  and  0.00)  these  grew  rapidly 
and  'occluded'  (several  free  boundary  drift  points  crossed  so  that  the 
free  boundary  formed  a  closed  loop)  shortly  after  six  inertial  periods. 
This  condition  caused  computational  difficulties  and  forced  termination 
of  the  runs. 
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Table  6 

Frontal  propagation  speeds  for  10  km  initial  disturbance 


Froude  Number 

Initial  Position 

-0.02 

5.5 

0.00 

7.5 

0.01 

7.5 

Speed  (cm/s) 
2.52 
-2.48 
-5.03 
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The  bulges  associated  with  the  initial  disturbance  moved  and 
amplified  over  this  time  span.  The  phase  speeds  and  amplitudes  were 
measured  by  the  same  method  that  was  discussed  in  chapter  5.  In  linear 
models  amplitude  growth  rates  have  an  exponential  time  dependence  which 
appears  as  an  imaginary  component  of  the  disturbance  phase  velocity 
(Orlanski,  1969).  We  assume  such  a  simple  relation: 

A  =  Aq  exp(gt)  (6.3.1 ) 

where  t  is  the  time,  AQ  is  the  initial  amplitude,  A  is  the  amplitude  at 
time  t  and  g  is  the  growth  rate.  The  formula  for  the  doubling  time 
derived  from  this  is: 


T2  *  in(2)/g  (6.3.2) 

where  T2  1s  the  doubling  time.  The  results  of  these  calculations  are 
shown  in  Table  7. 

The  effects  of  the  short  wavelength  disturbances  were  excluded 
from  the  amplitude  calculations  readily  by  a  curve  fitting  procedure. 

It  was  not  possible,  however,  to  entirely  exclude  these  effects  from 
the  phase  speed  calculations.  This  may  partially  explain  the  lack  of 
a  simple  pattern  in  the  phase  speed  results,  such  as  was  seen  for  the 
frontal  propagation  results  in  Tables  5  and  5.  The  results  are,  however, 
in  reasonable  agreement  with  the  observationally  determined  values  of 
4  to  8  cm/s.  The  doubling  times  agree  roughly  with  the  results  of  other 
models  (e.g.  Orlanski,  1969;  Lipps,  1963;  Tareev,  1965).  These  models 
were  linear,  small  amplitude  models.  The  observational  results  of 
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Table  7 

Results  for  10  km  initial  disturbance 


Froude  Number 

Phase  Speed 
(cm/s) 

Final  Amplitude 
(km) 

Doubling  Time 
(Inertial  Periods) 

-0.02 

7.18 

18.0 

7.08 

0.00 

11.32 

13.1 

15.40 

0.01 

3.60 

13.5 

13.86 
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Hansen  (1970)  indicate  that  these  times  are  much  too  short.  However, 
the  observational  results  were  based  primarily  on  larger  amplitude 
disturbances.  As  will  be  seen  in  section  6.4,  the  amplitude  of  the 
disturbance  can  have  a  profound  effect  on  the  growth  rate. 

The  negative  Froude  number  case  showed  the  greatest  tendency 
toward  instability  according  to  the  frontal  displacement  amplitude 
criterion  defined  in  chapter  5.  The  growth  of  the  main  disturbance  was 
much  faster  and  the  smaller  disturbances  began  earlier  and  were  more 
numerous.  The  Fr  =  0.01  case  was  the  least  unstable  in  this  respect. 

The  relative  stabilities  from  the  energy  calculations  may  be  used  for 
comparison.  The  mean  and  fluctuation  components  of  the  kinetic  and 
potential  energies  were  calculated  from  (2 . 6. 3)- (2 . S. 6) .  The  ratios  of 
fluctuation  to  mean  energies  are  shown  in  Table  8  for  the  initial  and 
final  (after  ten  inertial  periods)  values. 

These  results  showed  interesting  differences  from  those  in 
Table  7.  The  Fr  *  0.01  case  showed  the  greatest  increase  in  fluctuation 
kinetic  energy  and  the  Fr  =  -0.02  case  had  the  smallest  increase.  This 
would  imply  that  the  Fr  =  0.01  case  was  the  most  unstable.  The  po¬ 
tential  energy  ratios  also  show  a  different  pattern  from  the  frontal 
displacement  amplitude  changes,  with  F^  =  0.01  the  most  unstable  and 
Fp  -  -0.02  the  least.  The  F  =  0.01  case  had  a  loss  of  total  potential 
energy  while  the  other  two  increased.  The  Fp  =  -0.02  case  had  the 
smallest  loss  of  total  kinetic  energy  and  the  Fr  *  0.01  case  had  the 
greatest. 
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Table  8 

Fluctuation  and  mean  energy  for  10  km  initial 


Kinetic 

Initial 

Mean 

Fluctuation 
Ratio  (xl0~3) 
Final 
Mean 

Fluctuation 
Ratio  (xlO*3) 
Potential 
Initial 
Mean 

Fluctuation 
Ratio  (xlO-4) 
Final 


Froude  Number 

-0.02  0.00 


45.96  42.18 

0.2408  0.2415 

5.24  5.72 

37.32  33.39 

0.8927  1.1413 

23.6  34.2 


2393.  3007. 

0.4029  0.4091 

1.39  1.36 


Mean 

Fluctuation 
Ratio  (xlO*^) 


2940.  3026. 

0.3040  1.1019 

2.73  3.64 


disturbance 

0.01 

42.91 

0.2800 

6.53 

36.56 

2.5395 

69.5 

3201. 

0.4177 

1.31 

3188. 

2.0060 

6.29 
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For  small  amplitude  meanders  the  cross  stream  flow  has  a 
significant  effect  on  the  behavior  of  the  free  boundary  and  on  the 
energy  flow  in  the  outer  region.  Although  the  doubling  times  are  short 
compared  to  observational  results,  they  are  in  general  agreement  with 
other  small  amplitude  model  results.  The  negative  Froude  number  case 
appears  to  be  the  most  unstable  in  frontal  amplitude,  but  the  most 
stable  in  the  outer  region  flow.  It  also  matched  most  closely  the 
observed  phase  speed  and  had  the  most  rapid  amplitude  growth.  The  two 
stability  criteria  are  clearly  related  to  different  aspects  of  the  flow, 
since  the  results  for  the  two  are  opposite. 

6.4  -  Gulf  Stream  Large  Amplitude  Perturbations 

The  second  group  of  tests  which  were  run  for  the  Gulf  Stream 
front  had  a  large  amplitude  to  wavelength  ratio  with  a  value  of  0.333. 
The  initial  disturbance  for  these  cases  had  wavelengths  of  300  km,  as 
in  the  small  amplitude  perturbation  tests,  but  had  an  amplitude  of 
100  km.  The  same  set  of  cross  stream  Froude  numbers,  -0.02,  0.00,  and 
0.01  was  used.  All  of  the  cases  were  run  to  ten  inertial  periods  and 
the  Fr  =  -0.02  case  was  extended  to  thirty  four  inertial  periods,  as  a 
test  of  the  long  term  behavior  of  the  model.  An  additional  case  with 
Fr  *  -0.02  and  an  inverse  bulge  was  also  run.  The  inverse  bulge  was  a 
perturbation  which  extended  the  front  to  the  left  instead  of  to  the 
right  as  in  the  other  cases. 

The  frontal  propagation  speeds  were  calculated  in  the  manner 
described  in  section  6.3,  and  are  shown  in  Table  9.  The  inverse  bulge 
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Table  9 

Frontal  propagation  speed  for  100  km 

initial  disturbance 

Froude  Number 

Initial  Position 

Speed  (cm/s) 

-0.02 

3.5 

2.42 

-0.021 

12.5 

2.41 

0.00 

5.5 

-1.76 

0.01 

7.5 

-4.51 

. - 


case  has  a  Froude  number  designation  -0.021  to  distinguish  it  from  the 

'normal'  case  with  =  -0.02. 

r 

The  two  negative  Froude  number  cases  show  excellent  agreement 
with  each  other  and  with  the  results  in  Tables  5  and  6.  The  zero 
Froude  number  case  shows  a  significant  reduction  in  magnitude  from  the 
previously  observed  values.  The  Fr  »  0.01  speed  is  also  smaller  in 
magnitude  than  was  seen  earlier. 

Frontal  position  and  depth  contours  after  ten  inertial  periods 
are  shown  in  Fig.  12a-d).  Since  the  flow  is  primarily  geostrophic, 
the  movement  of  the  water  is  nearly  along  the  depth  contours.  These 
contours  therefore  approximate  streamlines.  The  contours  are  much 
closer  together  in  Fig.  12a)  than  in  b)  or  c).  This  indicates  that 
the  downstream  velocities  are  highest  for  Fr  *  -0.02.  The  fluctuations 
in  the  deepest  contours  near  the  bottom  of  each  figure  are  the  result 
of  waves  propagating  along  the  interface.  The  Fr  *  -0.02  case  appears 
to  be  much  less  affected  by  these  than  the  other  two  cases.  In  all  of 
the  figures  the  contours  are  closer  together  on  the  upstream  side  of 
the  main  bulge  (towards  the  bottom  of  the  figure)  than  on  the  downstream 
side.  This  indicates  increased  velocities  on  the  upstream  side.  The 
increased  momentum  seems  to  have  carried  the  water  farther  away  from 
the  unperturbed  frontal  position  as  it  flowed  around  the  main  bulge. 

The  points  of  maximum  deflection  are  progressively  farther  downstream 
for  the  deeper  contours.  The  deepest  contour  for  Fr  =  -0.02  appears  to 
have  the  greatest  angle  relative  to  the  unperturbed  frontal  position 
and  its  maximum  deflection  point  is  carried  further  downstream  than  in 
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the  other  two  cases.  In  all  three  cases  the  deepest  contour  showed 
approximately  the  same  amount  of  deflection  as  the  front.  Thus,  the 
cross  stream  flow  is  seen  to  affect  the  downstream  velocity  and  the 
direction  of  the  flow  across  the  entire  frontal  region.  The  inverse 
bulge  case  shown  in  Fig.  1 2d)  was  the  only  one  in  which  the  initial 
disturbance  did  not  amplify.  In  that  case  the  contours  show  the  most 
crowding  on  the  downstream  side  of  the  main  bulge.  The  deepest  contour 
showed  almost  no  deflection. 

Time  sequences  of  the  frontal  position  are  shown  in  Fig.  13a-d) 

for  Fr  =  -0.02,  0.00,  0.01,  and  -0.021.  The  sequences  start  with  the 

initial  curve  at  the  left.  Subsequent  frontal  positions  are  shown  at 

two,  six,  and  ten  inertial  periods.  All  of  these  cases  show  an  initial 

steepening  at  the  forward  side  of  the  bulge.  In  the  F^  =  -0.02  case  the 

steepening  abated  later  in  the  run  and  the  bulge  became  rounder.  In 

the  other  three  cases  the  steepening  continued,  particularly  in  the 

F  =  0.01  and  -0.021  cases, 
r 

In  each  case  a  complementary  bulge  (opposite  in  direction  from 
the  original  bulge)  formed  in  front  (downstream)  of  the  original  dis¬ 
turbance.  These  complementary  bulges,  in  turn,  had  complementary 
bulges  form  in  front  of  them.  Simmons  and  Haskins  (1979)  also  observed 
the  development  of  disturbances  downstream  of  an  initial  perturbation. 

In  their  experiment  the  downstream  disturbances  grew  to  a  size  equiv¬ 
alent  to  the  evolved  size  of  the  initial  perturbation.  The  movements 
of  the  complementary  bulges  suggest  group  velocities  much  greater  than 
the  phase  velocities  of  the  perturbations.  If  we  regard  these  distur- 


bances  as  dispersive  components  of  the  flow,  the  associated  group 
velocities  would  be  80  cm/s  for  =  -0.02,  66  cm/s  for  Fr  *  0.00,  and 
58  cm/s  for  F^  =  0.01.  The  complementary  bulges  developed  weakly  for 
F^  =  -0.021,  but  suggest  a  group  velocity  equal  to  the  Fr  =  -0.02  case. 
These  figures  agree  with  the  (widely  scattered)  group  velocities  in 
the  models  reviewed  by  Hansen  (1970)  which  ranged  from  11  to  100  cm/s. 

The  main  disturbance  wavelengths  in  all  cases  tended  to  grow. 

This  effect  was  most  pronounced  in  the  Fr  =  -0.02  test.  Two  explanations 
are  possible.  One  assumes  that  these  disturbances  are  dispersive,  and 
the  other  assumes  that  they  are  separable  into  spatially  and  temporally 
growing  components. 

The  initially  imposed  disturbance  has  tne  form  F(y)  =  B(0) 
[cos(ivy)  +  0.5],  where  8(0)  is  a  boxcar  function  (square  wave)  one 
wavelength  wide  and  centered  at  zero.  .  The  spatial  Fourier  transform 
of  this  is  F(k)  *  0.5[sinc(l-k)  +  sinc(l+k)  +  sinc(k)j,  where  k  is  the 
alongfront  wave  number  and  sinc(s)  =  sin(7rs)/irs.  The  initial  dis¬ 
turbance  is  thus  harmonically  rich.  If  these  waves  were  dispersive, 
then  as  the  various  components  begin  to  separate,  the  disturbed  region, 
and  hence  the  apparent  wavelength,  would  grow.  For  the  Froude  number 
to  affect  the  rate  of  increase  of  the  wavelength,  it  would  have  to 
affect  the  disoersive  character  of  the  perturbations. 

The  increase  in  wavelength  does  not,  however,  require  dispersive 
behavior.  As  suggested  by  Thacker  (1976),  the  disturbance  may  be 
separable  into  spatially  and  temporally  growing  components  (see  for 


example  his  eq.  18).  A  wave  which  has  both  spatially  and  temporally 
growing  components  would  show  an  increase  in  apparent  wavelength  as 
the  two  components  separated  in  space.  To  affect  the  rate  of  the 
apparent  increase  of  the  wavelength,  the  Froude  number  would  have  to 
affect  the  growth  rates  of  the  two  components  and/or  the  propagation 
rate  of  the  spatially  growing  component. 

The  shapes  of  the  bulges  in  Fig.  13  are  also  indicative  of  the 
nature  of  their  propagation.  The  trailing  sides  of  the  bulges  tended 
to  become  smoother  and  flatter  with  time.  Although  the  locations  of 
the  rearmost  parts  of  the  perturbations  became  difficult  to  distinguish 
precisely,  they  generally  retained  their  positions  or  advanced  slightly. 
The  forward  sides  and  the  amplifying  portions  of  the  perturbations 
clearly  moved  forward.  This  supports  the  hypothesis  that  the  distur¬ 
bances  are  both  spatially  and  temporally  growing. 

The  amplitudes,  phase  speeds  and  doubling  times  were  computed 
in  the  same  manner  as  in  section  6.3.  These  are  shown  in  Table  10. 

All  of  the  tests  showed  amplification  except  the  inverse  bulge. 
The  flow  in  the  outer  region  was  apparently  unable  to  continue  to  supply 
energy  to  this  perturbation.  The  disturbance  amplitudes  did  not  show 
a  temporally  uniform  growth  in  any  of  these  tests.  The  results  in 
Table  10  must  then  be  regarded  as  having  significant  possible  errors, 
so  that  there  may  be  little  real  difference  in  the  actual  doubling  times. 
This  would  imply  that  the  cross  stream  flow  does  not  have  a  major  effect 
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Table  10 

Results  for  100  km  initial  disturbance 


Froude  Number 

Phase  Speed 
(cm/s) 

Final  Amplitude 
(km) 

Doubling  Time 
(Inertial  Periods) 

-0.02 

3.94 

108.3 

86.9 

-0.021 

3.64 

93.0 

- 

0.00 

2.94 

107.2 

99.7 

0.01 


2.54 


110.0 


72.7 
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on  the  growth  rates  for  large  amplitude  perturbations,  while  it  does  for 
small  amplitude  disturbances. 

The  doubling  times  seen  here  are  much  longer  than  those  seen  for 

the  small  amplitude  disturbances  (the  amplification  rates  are  lower)  and 

are  much  closer  to  the  observationally  determined  values.  As  noted 

above,  the  observational  results  were  obtained  for  finite  amplitude 

meanders.  Thus,  the  tests  for  the  large  amplitude  disturbance  were  more 

nearly  representative  of  the  observational  conditions  than  were  the 

small  amplitude  tests.  These  results  indicate  that  the  nonlinear 

effects  in  finite  amplitude  disturbances  significantly  reduce  the  growth 

rate.  In  the  model  of  Orlanski  and  Cox  (1973),  the  growth  rate  for  the 

finite  amplitude  disturbances  was  an  order  of  magnitude  lower  than  that 

for  the  linear  (small  amplitude)  models.  The  small  amplitude  growth 

-6  -1  -6  -1 

rates  measured  here  were  0.5  x  10  s"  to  1.1  x  10"  s  versus  the 

0.96  x  10  s  to  1.2  x  10  s  rates  seen  by  Orlanski  and  Cox.  The 

-6  -1  -6 

finite  amplitude  disturbance  rates  were  0.081  x  10  s"  to  0.093  x  10 
s'1  here  «ersus  0.085  x  10"®  s"^  to  0.22  x  10“®  s"1,  showing  very  good 
agreement  between  the  two  models.  Simmons  and  Hoskins  (1978,  1979)  and 
Pedlosky  (1979)  also  found  a  severe  reduction  in  growth  rates  as  the 
instability  amplitude  increased. 

The  phase  speeds  were  also  lower  here  than  for  the  small  ampli¬ 
tude  disturbances.  These  speeds  are  slightly  lower  than  the  4  to  8  cm/s 
speeds  obtained  from  observations.  A  trend  of  increasing  speed  with 
decreasing  Froude  number  can  be  identified. 
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The  character  of  the  energy  changes  also  differed  significantly 
from  that  in  the  small  amplitude  tests.  The  fluctuation  to  mean  com¬ 
ponent  ratios  are  shown  in  Table  11. 

The  inverse  bulge  case  showed  a  decrease  in  both  the  fluctuation 
components,  which  is  consistent  with  the  reduction  in  the  amplitude  of 
the  disturbance  which  was  observed  in  Table  10.  The  other  three  tests 
showed  increases  in  the  ratios.  Although  not  directly  comparable,  the 
increases  were  smaller  than  those  seen  in  Table  8  for  the  small  ampli¬ 
tude  disturbances.  This  also  indicates  that  the  amplification  processes 
are  slower  for  the  large  amplitude  perturbations.  The  differences  seen 
between  the  three  Froude  number  tests  were  small.  This  agrees  with  the 
results  of  the  growth  rate  calculations  showing  that  the  effect  of  the 
Froude  number  was  greatly  reduced  for  these  tests. 

To  determine  the  long  term  behavior  of  the  model,  the  test 
with  F  =  -0.02  was  continued  for  a  total  time  of  thirty  four  inertial 
periods.  During  this  time  the  frontal  disturbance  continued  to  advance 
and  amplify.  Depth  contours  are  shown  in  Fig.  14  at  eighteen  and  thirty 
four  inertial  periods.  The  general  pattern  shown  in  Fig.  12  is  continued, 
with  the  wavelength  increasing,  the  original  bulge  becoming  less  steep 
and  the  complementary  bulge  in  front  of  it  becoming  deeper  and  more 
steep.  Due  to  the  periodic  boundary  conditions  and  the  relatively  high 
group  velocity,  the  disturbances  which  were  advancing  ahead  of  the  main 
perturbation  had  crossed  the  top  boundary,  reentering  the  grid  at  the 
bottom  boundary  and  had  caught  up  to  the  trailing  (bottom)  side  by  the 
end  of  thirty  four  inertial  periods.  At  this  point  they  were  beginning 
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Table  11 

Fluctuation  to  mean  energy  ratios  for  100  km 
initial  disturbance  amplitude 


-0.02 

Froude 

-0.021 

Number 

0.00 

0.01 

Kinetic  (xlO"3) 

Initial 

119. 

213. 

125. 

124. 

Final 

158. 

186. 

215. 

274. 

Potential  (xlO'4) 

Initial 

29.2 

60.9 

28.3 

28.9 

Final 

42.2 

22.1 

64.3 

70.2 
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to  affect  the  behavior  of  the  main  perturbation,  so  the  test  was 
terminated. 

The  propagation  of  the  front  across  the  grid  continued 
steadily,  with  the  frontal  propagation  speed  varying  less  than  5 %  from 
the  value  given  in  Table  10.  The  evolution  of  the  perturbation  also 
continued  at  approximately  the  same  rate  as  for  the  first  ten  inertial 
periods.  The  results  of  the  phase  speed  and  amplification  calculations 
are  shown  in  Table  12. 

As  noted  above,  the  increase  of  the  disturbance  amplitude  was 
not  uniform  in  time,  thus  creating  the  differences  in  the  doubling 
times.  These  results  are,  however,  reasonably  close  to  each  other 
and  still  show  a  marked  increase  from  the  small  amplitude  case.  As 
the  wavelength  of  the  main  bulge  increased,  the  kinetic  energy  fluctu¬ 
ation  to  mean  ratio  showed  non-uniform  behavior.  The  fluctuation  to 

_3 

mean  kinetic  energy  ratio  rose  from  153  x  10  at  ten  inertial  periods 

_3 

to  218  x  10  at  eighteen  inertial  periods.  It  then  declined  to 
178  x  10"3  at  thirty  four  inertial  periods.  This  may  have  been  a  result 
of  the  inertial  period  oscillations.  The  kinetic  energy  in  the  inertial 
period  oscillations  increased  with  time  as  shown  in  Fig.  15.  The 
oscillations  which  were  initially  small  began  increasing  significantly 
after  about  fourteen  inertial  periods.  Orlanski  and  Cox  (1973)  observed 
a  similar  increase  in  near  inertial  oscillations  in  their  model.  In 
that  study  the  inertial  period  waves  were  shown  to  be  caused  by  baro- 
clinic  instability. 


“>■* - 
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Table  12 

Results  for  extended  run  of  100  km  initial  disturbance 


Time 

Phase  Speed 

Final  Amplitude 

Doubling  Time 

(Inertial  Periods) 

(cm/s) 

(km) 

(Inertial  Periods) 

18 

2.64 

120.1 

68.12 

34 

2.91 

131 .0 

36.54 

Spectra  were  calculated  for  the  data  shown  in  Fig.  15  for  the 
time  intervals  0-14  and  14-34  inertial  periods,  and  are  shown  in  Fig. 

16.  The  error  bars  show  the  90"  confidence  intervals.  After  fourteen 
inertial  periods  the  energy  in  the  inertial  frequency  band  (and  its 
harmonics)  greatly  increases.  The  low  frequency  peak  in  16a)  is  missing 
entirely  in  16b).  In  general,  the  ea  iier  segment  showed  lower  variance, 
as  expected  from  Fig.  15.  The  slope  at  the  higher  frequencies  is 
approximately  -1/3.  Processes  involving  turbulence  usually  exhibit  a 
-5/3  slope.  This  behavior  is  consistent  with  the  model  operating  in  an 
inviscid  manner,  as  initially  intended. 

5.5  -  Sargasso  Sea  Short  Wavelength  Perturbation 

Persistent  oceanic  fronts  have  been  observed  in  the  Sargasso 
Sea  (Voorhis  and  Hersey,  1964;  Voorhis,  1969;  Katz,  1969).  A  set  of 
oceanographic  conditions  characteristic  of  this  region  were  used  to 
define  the  second  frontal  regime  considered  in  this  investigation.  The 
parameters  used  in  the  model  are  listed  in  Table  13. 

These  fronts  are  much  shallower  than  the  Gulf  Stream  front, 
and  are  also  narrower  (smaller  Rossby  radius).  The  maximum  internal 
wave  speed  is  lower  and  the  inertial  period  is  longer,  implying  that 
the  frontal  processes  will  be  slower. 

Observations  of  the  surface  thermal  field  (e.g.  Voorhis,  1969) 
indicate  the  presence  of  a  range  of  perturbation  wavelengths.  Since  the 
wavelength  of  a  disturbance  may  affect  its  stability  characteristics. 


Model  parameters  for  Sargasso  Sea  front 


Stratification  (ap/p  ) 

oo 

Latitude 

°b 

Do 

Inertial  Period 
g‘  (reduced  gravity) 
c  (maximum  internal  wave  speed) 

Maximum  internal  wave  speed  at  the  front 
X  (baroclinic  Rossby  radius) 


0.CO13 
30°  N 
50  m 
200  m 
23.94  hr 
0.0127  m/s 
1.596  m/s 
0.798  m/s 
21.9  km 
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tests  were  performed  with  two  initial  disturbance  wavelengths.  These 
were  chosen  as  30  km  and  100  km  to  cover  the  observed  range.  All  of 
the  tests  were  run  for  ten  inertial  periods. 

For  the  short  wavelength  tests  a  grid  spacing  of  2  km  and  a 
time  step  of  359.129  s  were  used.  This  results  in  a  Courant  number  of 
0.81,  which  is  very  close  to  the  0.82  value  used  previously.  A  small 
initial  amplitude  of  3  km  was  used.  To  determine  the  dependence  of 
the  frontal  evolution  on  the  cross  stream  flow,  the  same  set  of  Froude 
numbers  that  was  used  in  the  previous  tests  was  also  used  here.  The 
frontal  propagation  speeds  were  calculated  as  before  and  are  shown  in 
Table  14. 

As  seen  previously,  the  Froude  number  is  a  major  determinant 

of  the  propagation  speed.  The  speeds  seen  here  are  much  smaller  than 

those  for  the  Gulf  Stream  front,  consistent  with  the  expectation  of 

generally  slower  movement.  The  speed  variation  was  approximately 

0.5  cm/s  for  a  0.01  increment  in  F„. 

r 

Time  sequences  of  the  frontal  position  at  one  inertial  period 
intervals  are  shown  in  Fig.  17a-c)  for  Fr  *  -0.02,  0.00,  and  0.01, 
respectively.  The  x  coordinate  has  been  expanded  by  a  factor  of  five 
in  all  cases.  The  evolution  of  these  disturbances  was  initially  similar 
to  the  evolution  of  the  Gulf  Stream  front  disturbances.  An  initial 
steepening  of  the  leading  edge  and  the  formation  of  a  complementary 
bulge  occurred  for  all  cases.  However,  for  the  Fr  =  0.00  and  0.01  cases 
the  comolementary  bulge  rapidly  diminished  in  amplitude  and  was  an 
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indistinct  feature  within  five  inertial  periods.  Because  of  this  rapid 
dimunition,  meaningful  group  velocities  could  not  be  calculated.  For 
the  =  -0.02  case  the  complementary  bulge  was  maintained  longer,  but 
eventually  vanished  also.  Based  on  the  motion  during  the  first  eight 
inertial  periods  the  group  velocity  was  8.08  cm/s  for  this  case.  This 
is  much  lower  than  the  group  velocities  seen  in  the  Gulf  Stream  large 
amplitude  tests.  The  F^  =  0.01  case  developed  short  wavelength  dis¬ 
turbances,  similar  to  those  seen  in  the  Gulf  Stream  small  amplitude 
tests,  although  here  they  did  not  grow  to  occlusion.  The  phase  speeds 
and  final  amplitudes  are  shown  in  Table  15. 

The  two  non-zero  Froude  number  cases  had  similar  phase  speeds, 
while  the  phase  speed  for  the  F^  =  0.00  case  was  somewhat  lower.  Al¬ 
though  there  are  no  reliable  observational  data  for  comparison,  these 
results  seem  to  be  in  a  reasonable  range.  They  are  generally  smaller 
than  the  speeds  for  the  Gulf  Stream  small  amplitude  tests.  The  Froude 
number  does  not  seem  to  have  a  major  effect  on  these  speeds. 

All  of  the  disturbances  showed  a  decrease  in  amplitude.  Con¬ 
sidering  that  the  losses  of  total  energy  were  greater  here  than  for  the 
Gulf  Stream  front,  it  appears  that  these  runs  may  have  been  overdamped 
by  interfacial  friction.  The  same  value  of  Cg  was  used  here  as  for 
the  Gulf  Stream  front.  However,  since  this  flow  was  shallower  and 
slower,  a  lower  value  may  have  been  more  appropriate.  In  spite  of  this, 
a  trend  in  the  results  is  clear.  The  Fr  =  -0.02  case  maintained  a 
disturbance  closest  to  the  original  amplitude.  The  amplitudes  of  the 
other  two  tests  decreased  steadily  to  about  60“  of  the  initial  value. 
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Table  15 

Results  for  30  km  wavelength  initial  disturbance 


Phase  Speed 


Froude  Number  (cm/s) 

-0.02  4.69 

0.00  3.62 

0.01  4.59 


Final  Amplitude 
(km) 

1.98 

1.20 

1.12 


This  indicates  the  greatest  tendency  towards  stability  in  the  positive 
Froude  number  case,  as  was  seen  above,  when  using  the  amplitude 
stability  criterion. 

The  changes  in  the  fluctuation  to  mean  energy  ratios  were 
also  similar  to  those  of  the  previous  tests.  These  are  shown  in 
Table  16. 

These  results  again  indicate  that  the  negative  Frouae  number 
case  has  the  most  stable  flow  in  the  outer  region  and  the  most  unstable 
disturbance  amplitude. 

6.6  -  Sargasso  Sea  Long  Wavelength  Perturbation 

As  mentioned  above,  a  set  of  tests  were  run  for  the  Sargasso 
Sea  front  with  an  initial  perturbation  wavelength  of  100  km.  For  these 
tests  the  grid  spacing  was  4  km  and  the  time  step  was  71 S. 258  s,  which 
again  results  in  a  Courant  nun®er  of  0.81.  Since  the  Froude  number 
dependence  had  been  closely  investigated  previously,  a  single  value 
of  0.00  was  used  for  these  tests.  Two  model  runs  were  made  with  initial 
amplitudes  of  3  and  30  km.  The  amplitude  to  wavelength  ratios  were  then 
0.03  and  0.30,  respecti vely.  The  frontal  propagation  speeds  are  shown 
in  Table  17. 

These  speeds  are  significantly  higher  than  those  for  the  short 
wavelength  case.  This  and  the  large  difference  between  these  two  cases 
indicate  a  sensitivity  to  the  wavelength  of  the  perturbation. 
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Table  16 

Fluctuation  to  mean  energy  ratios  for  30  km  wavelength 
initial  disturbance 


Kinetic  (xlO"^) 
Initial 
Final 

Potential  (x!0“4) 
Initial 
Final 


Froude  Number 

-0.02 


.404 

1.16 

.170 

.586 


0.00 

.319 

8.17 

.178 

3.23 


0.01 

.401 

9.78 

.165 

6.09 
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Table  17 

Frontal  propagation  speeds  for  100  km  wavelength 
initial  disturbance 


Initial  Amplitude 

Initial 

Speed 

(km) 

Position 

(cm/s) 

3 

7.5 

-1.95 

30 

7.5 

-1.51 

-  » 
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A  time  sequence  of  the  frontal  position  for  the  small  amplitude 
case  is  shown  in  Fig.  17.  The  x  coordinate  has  been  expanded  by  a  factor 
of  five.  The  time  sequence  for  the  large  amplitude  case  is  shown  in 
Fig.  19.  The  x  coordinate  is  not  expanded  in  that  figure. 

The  leading  edge  of  the  small  amplitude  disturbance  followed 
the  general  pattern  of  initial  steepening,  formation  of  a  complementary 
bulge,  and  then  a  period  of  decreasing  steepness  at  the  leading  edge 
and  a  deepening  of  the  complementary  bulge.  In  that  respect  this 
perturbation  benaved  more  like  the  Gulf  Stream  cases  than  the  Sargasso 
Sea  short  wavelength  case  where  the  complementary  bulge  did  not  remain 
a  permanent  feature.  The  trailing  edge  was  anomalously  active, 
initially  advancing,  then  retreating,  and  then  advancing  again.  During 
the  retreating  phases  it  moved  upstream,  became  very  steep,  and  even¬ 
tually  collapsed.  This  type  of  motion  was  not  seen  along  the  Gulf 
Stream  front  in  sections  6.3  and  6.4.  If  this  process  actually  occurs 
in  the  Sargasso  Sea,  the  meander  evolution  would  be  more  complex,  making 
accurate  observation  of  phase  speeds  and  growth  rates  more  difficult. 

The  large  amplitude  perturbation  (Fig.  19)  followed  the  general 
pattern  seen  in  the  Gulf  Stream  cases.  This  test  also  appeared  to  show 
the  overdamping  mentioned  above.  In  addition  to  the  lack  of  disturbance 
amplification  there  was  significantly  less  activity  along  the  interface. 
Fig.  20  shows  interfacial  depth  contours  at  the  initial  and  final  stages. 
Although  the  general  pattern  is  similar  to  that  in  Fig.  12b),  the  very 
short  wavelength  fluctuations  along  the  deeper  contours  are  much  smaller 
in  Fig.  20b). 
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The  group  and  phase  speeds  and  the  doubling  time  results  are 
shown  in  Table  18. 

The  phase  speeds  were  lower  than  for  the  short  wavelength  case 
and  the  group  speeds  were  higher  than  the  value  found  for  Fr  =  -0.02. 
Although  the  effects  of  the  possible  overdamping  are  not  clear,  these 
tests  clearly  show  greater  disturbance  amplification  tendencies.  The 
short  wavelength,  zero  Froude  number  test  showed  significant  dimunition, 
whereas  these  tests  showed  neutral  stability  for  the  30  km  case  and 
amplification  for  the  3  km  case.  The  3  km  case  had  a  doubling  time 
comparable  to  the  Gulf  Stream  small  amplitude  case. 

The  phase  and  group  speeds  showed  that  the  small  amplitude 
case  had  more  rapid  movement  than  the  large  amplitude  case,  as  was 
seen  for  the  Gulf  Stream.  The  amplification  rates  agreed  with  the 
prior  results  in  shewing  greater  amplitude  instability.  The  fluctuation 
to  mean  energy  ratios  shown  in  Table  19  also  reflect  that  pattern. 

The  growth  of  fluctuation  kinetic  energy  was  slower  for  the 
large  amplitude  test,  as  was  seen  before  in  the  Gulf  Stream.  The 
decrease  in  the  potential  energy  ratio  for  the  large  amplitude  case 
suggests  a  straightening  of  the  interfacial  depth  contours,  which  is 
supported  by  Fig.  19.  Comparing  the  3  km  amplitude  case  here  with  the 
short  wavelength  test  discussed  above  shows  that  the  growth  in  the 
kinetic  energy  ratio  for  the  long  wavelength  case  was  much  smaller 
than  that  for  the  short  wavelength  case  while  the  disturbance  amplitude 
growth  was  larger.  This  is  the  same  pattern  as  was  seen  above. 
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Table  18 

Results  for  100  km  wavelength  initial  disturbance 


Phase  speed  (cm/s) 

Group  speed  (cm/s) 

Final  amplitude  (km) 

Doubling  time 

(Inertial  periods) 


Initial  Amplitude  (km) 
3  30 

2.27  1.89 

17.4  14.6 

4.60  29.7 

6.95 


F/e  8/10 


A0-A115  519  DELAWARE  UNIV  NEWARK  COLL  OF  MARINE  STUOIES 

OCEANIC  FRONTAL  STABILITY*  A  NUMERICAL  MODEL’ <U) 

NOV  81  W  W  MARTIN  N00014-75-C-071* 

UNCLASSIFIED  C MS-0 2-81  NL 
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Table  19 

Mean  to  fluctuation  energy  ratios  for  100  km  wavelength 
initial  disturbance 


Initial 

Amplitude  (km) 

3 

30 

Kinetic  (xlQ~^) 

- 

Initial 

1.35 

54.3 

Final 

9.33 

63.4 

Potential  (xlO-4) 

Initial 

.646 

29.0 

Final 


1.35 


23.7 


CONCLUSIONS 


A  stable,  convergent  model  has  been  developed  to  study  the 
evolution  of  mesoscale  meanders  of  oceanic  fronts.  This  model  has 
focused  on  some  of  the  parameters  which  affect  the  stability  charac¬ 
teristics  of  these  meanders.  Oceanographic  conditions  representative  of 
two  regions,  the  Gulf  Stream  downstream  of  Cape  Hatteras  and  the  Sargasso 
Sea  at  30°  N,  were  used  to  set  several  of  the  basic  model  parameters. 

The  results  of  the  model  were  verified  where  possible  by 
comparison  with  observations.  The  only  observations  available  for 
comparison  with  the  model  results  are  from  large  amplitude  meanders 
of  the  Gulf  Stream.  Both  the  phase  speeds  and  amplitude  doubling 
times  computed  from  the  model  results  appropriate  for  these  tests 
were  in  satisfactory  agreement  with  the  observed  values. 

Since  the  model  employed  nearly  inviscid  dynamics,  it  should 
have  yielded  flows  which  nearly  conserved  potential  vorticity.  Stommel 
(1976)  gave  a  simple  model  for  the  Gulf  Stream  which  assumed  geostrophic 
balance  for  the  alongfront  velocity  and  uniform  potential  vorticity 
across  the  stream.  The  initial  conditions  of  the  present  model  also 
used  these  assumptions,  so  that  the  potential  vorticity  was  initially 
uniform  across  the  model  grid.  For  the  tests  which  ran  for  ten  inertial 
periods,  the  potential  vorticity  field  remained  uniform  to  within  3". 
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For  the  extended  test,  which  was  run  to  thirty  four  inertial  periods, 
the  potential  vorticity  remained  uniform  to  within  7%.  This  indicates 
that  the  model  was  conserving  potential  vorticity  to  within  acceptable 
limits.  The  small  discrepancies  are  easily  accounted  for  by  the  slight 
interfacial  friction  used  for  stability.  The  model  potential  vorticity 
field  thus  provides  an  independent  check  on  the  model  results. 

Several  scalings  were  used  in  the  development  of  this  model. 
These  were:  the  baroclinic  Rossby  radius  for  length,  the  inertial 
period  for  time,  and  the  maximum  internal  wave  speed  for  velocity. 

The  model  results  indicate  that  the  length  scale  was  chosen  properly, 
but  that  the  time  and  velocity  scales  were  not.  The  characteristic 
time  scale  for  the  meanders  was  on  the  order  of  many  inertial  periods, 
so  that  the  time  scaling  used  was  too  short  and  the  velocity  scaling 
too  large.  The  velocity  scaling  was,  however,  proper  for  the  fluid 
velocity  itself.  It  is  still  unclear  what  the  meander  time  scale 
should  be  in  terms  of  flow  parameters.  Although  this  model  and  others 
with  similar  scalings  yield  reasonable  results,  the  lack  of  a  proper 
time  scale  indicates  that  some  of  the  essential  physics  have  yet  to  be 
understood.  This  problem  will  most  likely  require  further  development 
of  theoretical  models  for  its  resolution. 

Tests  were  conducted  with  this  model  to  determine  the  effects 
of  variation  of  several  model  parameters  on  the  evolution  of  mesoscale 
meanders.  Of  those  which  can  affect  the  behavior  of  these  disturbances. 


the  wavelength/Rossby  radius  ratio,  amplitude/wavelength  ratio,  and 
cross  stream  Froude  number  were  investigated. 

One  of  the  main  purposes  of  the  many  linear  models  has  been  to 
deduce  the  wavelength  of  the  fastest  growing  instability.  Nevertheless, 
in  this  model,  since  the  predominant  wavelength  is  generally  known  from 
observation,  it  was  treated  as  a  given  model  parameter.  In  the  simple 
experiment  that  was  performed  to  investigate  the  effect  of  two  different 
wavelengths  on  stability,  significant  differences  were  seen.  While  the 
small  wavelength/Rossby  radius  disturbance  showed  a  decrease  in  ampli¬ 
tude,  the  large  wavelength/Rossby  radius  ratio  disturbance  showed 
significant  growth.  Furthermore,  a  large  wavelength/Rossby  radius  ratio 
disturbance  together  with  a  larger  wavelength/ampl itude  ratio  also 
showed  an  amplitude  decrease,  but  the  decrease  was  not  as  large  as  for 
the  short  wavelength  case.  For  the  Sargasso  Sea  conditions  a  30  km  long 
disturbance  was  seen  to  be  more  stable  than  a  100  km  long  disturbance. 
These  results  indicate  that  wavelength/Rossby  radius  ratio  is  a  signif¬ 
icant  parameter  for  this  model  which  affects  stability. 

The  initial  amplitude/wavelength  ratio  was  a  second  parameter 
shown  to  have  a  strong  influence  on  the  amplitude  growth  rate.  Several 
tests  were  made  with  amplitude  to  wavelength  ratios  from  0.03  to  0.33, 
using  both  sets  of  oceanographic  conditions.  The  small  amplitude  cases 
all  showed  more  rapid  amplitude  growth  than  the  large  amplitude  cases. 

For  the  small  amplitude  tests  using  Gulf  Stream  conditions,  the  growth 


rates  were  similar  to  those  given  by  several  prior  small  amplitude 
analytic  models  (Lipps,  1963;  Tareev,  1965;  Orlanski,  1969).  The  growth 
rates  obtained  from  the  corresponding  large  amplitude  tests  closely 
matched  the  rates  obtained  from  field  and  satellite  observations,  which 
were  primarily  of  large  amplitude  disturbances  along  the  Gulf  Stream 
(Hansen,  1970;  Halliwell  and  Mooers,  1979).  The  amplitude  growth  rates 
for  the  large  amplitude  disturbances  were  typically  an  order  of  magni¬ 
tude  smaller  than  those  for  the  small  amplitude  cases.  The  non-linear 
terms  in  the  governing  equations  likely  become  a  major  factor  in  the 
finite  amplitude  dynamics,  reducing  the  amplitude  instability  as  the 
amplitude  increases. 

The  cross  stream  Froude  number  was  shown  to  be  a  third  important 
parameter  for  the  stability  characteristics.  The  movement  of  the  front 
in  the  cross  stream  direction  was  related  linearly  to  this  parameter. 

A  negative  Froude  number  flow,  resulting  from  downward  entrainment  in 
the  inner  region,  was  required  to  eliminate  cross  stream  frontal  mi¬ 
gration  over  a  long  period  of  time.  This  downward  entrainment  also 
resulted  in  a  decrease  in  the  amplitude  doubling  time.  This  decrease 
was  by  as  much  as  a  factor  of  two,  compared  with  a  case  with  no  cross 
flow,  in  the  Gulf  Stream  small  amplitude  tests. 

Two  stability  criteria  were  used  in  these  tests.  One  was  based 
on  the  amplitude  of  the  disturbance  and  the  other  on  the  fluctuation  to 
mean  component  ratios  for  the  kinetic  and  potential  energies  of  the  flow 
in  the  outer  region.  These  two  criteria  were  shown  to  be  in  opposition. 


In  a  given  set  of  tests,  the  case  which  was  most  stable  by  one  crite¬ 
rion  was  the  least  stable  by  the  other.  This  suggests  that  the  flow 
in  the  outer  region  and  the  deformation  of  the  front  were  competing 
for  energy  from  the  same  source.  The  energy  source  most  likely  was 
the  potential  energy  in  the  parent  pool,  although  the  energy  analysis 
performed  was  not  adequate  to  determine  this  clearly.  The  transfer  of 
energy  was  affected  by  the  Froude  number,  wavelength,  and  particularly 
by  the  amplitude  of  the  disturbance,  although  these  effects  could  not 
be  quantified.  An  important  extension  of  this  model  would  be  a  proper 
determination  of  the  rates  and  pathways  of  the  energy  transfers. 

Inertial  period  oscillations  were  found  to  be  a  major  element 
in  the  temporal  variation  of  the  kinetic  and  potential  energies  of  the 
system,  although  they  did  not  appear  to  seriously  affect  the  frontal 
evolution  at  longer  than  inertial  periods.  These  oscillations  grew 
particularly  large  in  the  latter  portion  of  the  extended  Gulf  Stream 
test.  Orlanski  and  Cox  (1977)  also  found  large  growth  of  inertial 
period  oscillations  near  the  end  of  their  model  runs.  Their  results 
and  the  present  ones  suggest  that  oceanic  fronts  are  sensitive  to 
excitation  at  this  frequency. 

The  inertial  motions  were  particularly  intense  where  the  water 
velocity  was  the  highest,  near  the  inner/outer  boundary.  The  inner 
region  is  a  complex  zone  of  turbulent  mixing  and  interleaving  and  may 
be  able  to  dissipate  the  (relatively)  high  frequency  energy  in  the 
inertial  oscillations  and  prevent  their  growth  to  large  amplitude.  In 
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this  investigation  the  inner  region  was  modeled  very  simply,  with  no 
turbulent  energy  dissipation  occurring  within  it.  Kinetic  energy  was 
simply  advected  across  the  inner/outer  boundary  for  the  non-zero  Froude 
number  cases,  but  the  transfer  rates  were  controlled  by  the  inviscid 
flow  in  the  outer  region,  not  by  the  inner  region  flow.  The  effect  of 
this  simple  treatment  would  likely  be  most  important  for  positive  Froude 
number  flows,  when  water  is  flowing  from  the  inner  region  toward  the 
outer  region.  The  inclusion  of  an  inner  region  with  its  own  dynamics 
would  be  an  important  extension  of  this  model. 

Whatever  the  mechanism  for  its  dissipation,  energy  at  inertial 
frequencies  remains  important.  Future  modeling  efforts  should  be 
designed  with  the  realization  that,  unless  dispersed,  inertial  motions 
can  grow  to  dominate  the  flow.  Corresponding  field  observations  to 
determine  the  rate  of  energy  dissipation  in  the  inner  region  and  the 
direction  of  the  flow  across  the  inner/outer  boundary  (the  sign  of  the 
Froude  number)  are  also  needed  for  a  better  understanding  of  the  dy¬ 
namics  of  this  region. 

These  results  show  the  utility  of  the  division  of  the  frontal 
zone  into  two  regions,  turbulent  inner  and  inviscid  outer,  as  suggested 
by  Garvine  (1979).  The  characteristics  of  the  flow  in  the  outer  region, 
such  as  the  conservation  of  potential  vorticity,  are  the  result  of 
inviscid  dynamics.  However,  without  a  turbulent  inner  region  the 
entrainment  (upward  or  downward)  necessary  for  a  cross  stream  flow 
cannot  occur.  This  region  is  also  potentially  needed  as  an  energy 
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sink  for  inertial  period  oscillations.  Extensions  of  the  present  work 
to  include  a  complete  energy  analysis  and  an  account  of  the  inner  region 
dynamics  would  be  useful  in  determining  the  overall  importance  of  this 
region.  Field  studies  are  also  needed  to  determine  the  direction  and 
magnitude  of  the  entrainment  and  the  size  of  the  turbulent  inner  region. 


APPENDIX  I  -  COMPUTATIONAL  IMPLEMENTATION 


Al.l  -  Finite  Differencing 

The  numerical  algorithm  used  here  is  a  finite  difference  scheme 
based  on  that  proposed  by  Lax  and  Wendroff  (1960)  for  the  solution  of 
hyperbolic  partial  differential  equations  in  conservation  form.  The 
essential  idea  of  this  method  is  the  expansion  of  the  initial  equations 
in  a  Taylor  series,  usually  with  second  order  accuracy. 

Since  strong  spatial  gradients  are  present  in  the  initial  states 
of  all  of  the  cases  dealt  with  here,  difficulties  were  encountered  with 
the  numerical  scheme.  An  unmodified  Lax-Wendroff  approach  becomes 
unstable  due  to  the  growth  of  the  second  order  terms  in  the  regions  of 
strong  gradients.  Similar  difficulties  were  encountered  with  several 
other  schemes  which  were  tested  (some  of  which  were  variants  of  Lax- 
Uendroff)  because  they  were  not  designed  to  deal  with  the  strong  gradients. 
The  method  finally  chosen  is  based  on  an  algorithm  by  Harten  and  Zwas 
(1972)  intended  for  computation  of  compressible  flows  which  include 
shocks.  It  is  a  Lax-Wendroff  variant  which  is  second  order  in  smoothly 
varying  regions  and  is  reduced  to  first  order  in  shock  regions.  The 
change  from  first  to  second  order  accuracy  is  controlled  by  a  variable, 

3,  which  Harten  and  Zwas  call  an  'automatic  switch1.  This  variable  is 
of  order  one  near  a  shock  and  approaches  zero  in  smooth  regions.  The 


second  order  terms  in  a  Lax-Wendroff  scheme  are  multiplied  by  the  comple¬ 
ment  of  the  automatic  switch  (1  -  0).  Thus  in  smooth  regions  where  3  is 
small,  the  scheme  will  be  second  order.  Near  shocks,  where  0  approaches 
one,  the  second  order  terms  will  be  diminished  or  eliminated  so  that  the 
scheme  becomes  first  order.  By  reducing  the  relative  size  of  the  second 
order  terms,  some  accuracy  is  lost  near  the  shock.  However,  the  solutions 
remain  convergent  and  are  not  subject  to  the  instabilities  introduced  at 
the  shock  by  the  second  order  terms  used  in  other  methods.  The  fact 
that  the  present  work  deals  with  fronts,  which  have  mild  spatial  gradients 
when  compared  to  shocks,  does  not  affect  the  applicability  of  this 
method.  However,  some  problems  are  encountered  with  the  form  of  the 
smoothing  or  'artificial  viscosity'  originally  suggested  by  Harten  and 
Zwas,  which  shall  be  discussed  below. 

The  'automatic  switch*,  0,  has  a  value  near  zero  in  regions  of 
small  spatial  gradients,  and  rises  to  a  maximum  of  unity  in  the  frontal 
region  near  the  inner/outer  boundary.  A  derivation  is  given  in  Harten 
and  Zwas  (1972).  The  form  used  in  this  model  for  the  x  and  y  direction 
switches  is: 


*  "r2J»  J 


,-2 


0y  ,  =  0.5[ j 0  ,  -  <5 ■  , |/max|o -  ,  -  5.  ,jj 

T,j4  i.j4 


(Al.lb) 


where  max  j  j  is  the  maximum  difference  in  the  x  (or  y)  direction.  The 


subscripts  i  and  j  refer  to  the  x  and  y  axis  coordinates  of  the  grid 
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points.  The  subscript  i  denotes  the  entire  range  of  the  grid  along  the 
x  or  y  axis. 

The  partial  differential  equation  to  be  solved,  written  in 
conservation  form  is,  from  (2.4.1) 

wt  =  fx  +  gy  +  2 

where  w  is  a  column  vector  and  f,  g,  and  z  are  functions  of  w.  Employing 

the  Jacobi ans  A  *  and  8  =  we  have 

3w  3w  , 
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The  terms  on  the  right  hand  side  of  (AT .2)  are  the  first  order 

terms  and  those  on  the  right  hand  side  of  (AT. 5)  are  the  second  order 

terms  of  the  Lax-Wendroff  scheme.  Representing  f.  ,  .  as  F+,  f.  .  .  as 

F',  as  G+,  etc. ,  and  using  FC+  as  +  f,+1J  -  - 

f.  ,  ./ 4,  etc.,  the  scheme  is 
'  *  •  *  J 


=  w 


1  ej 


At 

2h 


(F+  -  F')  +  1^-  (G+  -  G')  +  Atz 


At 

2h 


1  ('  - 
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C[^j-  (F+-F  +G+-G)+z] 


(A1.6) 


where  At  is  the  time  step  and  h  is  the  grid  spacing  which  is  the  same 
along  the  x  and  y  axes. 
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A1.2  -  Smoothing 

Virtually  all  finite  difference  algorithms  for  fluid  dynamic 
equations  employ  some  form  of  smoothing  or  'artificial  viscosity'.  The 
purpose  of  the  smoothing  is  more  a  matter  of  computational  stability 
than  an  attempt  to  model  actual  viscous  effects.  If  fluid  viscosity  is 
to  be  included  in  a  model  it  is  usually  incorporated  as  an  additional 
term  in  the  difference  equations.  Smoothing  effects  are  usually  intro¬ 
duced  by  comparing  the  calculated  value  of  a  variable  at  a  grid  point  to 
an  idealized  value.  The  idealized  value  is  commonly  taken  to  be  the 
average  of  the  four  nearest  neighbors.  A  smoothing  term  is  computed  as 
a  proportion  of  the  deviation  of  the  calculated  value  from  the  idealized 
value.  The  smoothed  value,  which  is  the  final  result,  is  the  sum  of  the 
calculated  value  and  the  smoothing  term.  In  most  algorithms  the  proportion 
of  the  deviation  of  the  calculated  value  from  the  idealized  value  used 
to  form  the  smoothing  term  is  the  same  for  all  grid  points.  The  Harten 
and  Zwas  scheme  uses  the  automatic  switch  to  vary  this  proportion  so 
that  the  smoothing  is  strongest  near  the  shock.  The  idealized  value  is 
still  obtained  by  averaging.  Although  this  method  is  acceptable  for 
shock  computations  and  for  most  front  models,  it  could  not  be  used  in 
the  present  study  due  to  incompatibility  of  the  linear  averaging  with 
the  curvature  of  the  interfacial  depth  and  velocity  fields. 

In  shock  computations  the  important  results  are  the  location  and 
velocity  of  the  shock  at  a  given  time.  The  details  of  the  computed  flow 
within  the  actual  shock  region  are  not  considered  important,  as  the 
shock  is  expected  to  form  on  a  small  scale  in  the  cross-shock  direction. 


Ideally,  the  shock  should  be  sub-grid  scale  to  represent  a  discontinuity. 
Thus,  in  shock  regions  the  error  caused  by  any  inaccuracies  in  the 
smoothing  is  not  of  great  importance.  Outside  of  the  shock  region  the 
flow  is  expected  to  be  either  uniform  or  slowly  varying.  In  these 
regions  the  magnitude  of  the  smoothing  term  is  small.  Therefore,  models 
dealing  with  shock  flows  which  usually  employ  this  averaging  technique 
in  the  smoothing  are  not  presented  with  significant  problems. 

Numerical  models  of  fronts  are  also  able  to  use  this  method  of 
smoothing  without  difficulty  due  to  the  assumption  of  a  linear  depth 
profile  which  is  usually  made.  In  this  assumption  the  depth  (or  height) 
of  the  interface  is  taken  to  increase  linearly  with  ai stance  away  from 
the  surface  expression  of  the  front.  If  an  idealized  value  is  computed 
by  averaging,  it  will  lie  along  the  linear  profile.  With  the  idealized 
value  being  equal  (or  nearly  so)  to  the  initially  assumed  value,  the 
smoothing  terms  will  act  as  a  force  which  restores  the  smoothed  value 
towards  the  assumed  profile  (Fig.  21a). 

In  the  present  case  the  interface  profile  has  an  exponential 
variation.  If  an  idealized  value  for  the  depth  at  a  given  grid  point  is 
calculated  by  averaging,  it  will  not  agree  with  the  initially  assumed 
value,  but  will  be  towards  the  'inside'  of  the  profile  curve.  A  smoothing 
term  based  on  this  type  of  idealized  value  will  bias  the  smoothed  value 
away  from  the  initially  assumed  exponential  curve  and  towards  a  straight 
line  (Fig.  21b).  Unlike  the  shock  computations,  the  details  of  the  flow 
in  the  rapidly  changing  region  are  of  great  importance,  and  unlike 


previous  frontal  models  the  averaging  technique  for  obtaining  idealized 
values  does  not  work.  Thus,  we  must  derive  a  method  for  calculating 
idealized  values  which  will  fit  the  assumed  curve  of  the  depth  profile 
more  accurately  and  will  not  introduce  bias.  The  exponential  curve 
(4.1.10)  used  for  the  initial  depth  values  is  a  useful  shape  for  use  in 
smoothing. 

$(ns)  =  1  -  (1  -  $b)e  s 

At  each  time  step  this  equation  is  used  to  predict  the  depths 
along  each  x-directed  grid  line  from  the  current  local  front  position. 
These  idealized  values  are  subtracted  from  the  calculated  depths  to  form 
an  array  of  deviations.  The  smoothing  term  at  each  grid  point  is 
proportional  to  the  difference  between  the  deviation  at  that  point  and 
the  average  of  the  deviations  at  the  two  neighboring  points. 

This  metnod  provides  a  consistent  smoothing  without  introducing 
a  bias.  The  basic  technique  may  also  be  used  for  the  along  front  compo¬ 
nent  of  velocity.  The  velocity  is  related  to  the  depth  by  (4.2.5) 

v  -  1  -  5 

The  same  procedure  of  smoothing  according  to  deviations  from  a 
predicted  profile  is  applied  along  x  directed  grid  lines  for  the  along 
front  velocity  to  provide  a  non-biasing  smoothing  operator.  Along  y 
directed  grid  lines  the  curves  are  expected  to  be  sufficiently  flat  so 
that  linear  smootnino  may  be  used  for  both  depth  and  along  frontal 
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The  component  of  velocity  in  the  across  front  direction  will 
have  small  values  and  small  deviations  from  a  linear  profile  so  that  the 
technique  based  on  averaging  neighbors  will  produce  small  smoothing 
terms.  This  component  of  velocity  will  also  be  very  sensitive  to  small 
changes  in  the  depth  and  the  other  velocity  component  so  that  a  smoothing 
method  more  sophisticated  cannot  be  justified.  Thus,  the  method  of 
averaging  nearest  neighbors  to  obtain  idealized  values  is  used  for  this 
variable. 

A1.3  -  Extrapolation  to  Inner  Region  Points 

The  scheme  detailed  in  (A1.6)  requires  values  for  the  three 
variables  at  four  neighboring  points.  This  presents  no  difficulty  in 
the  center  of  the  grid  or  at  the  three  boundaries  discussed  in  section 
2.4.  Points  near  the  free  boundary  may  have  nearest  neighbors  which  are 
in  the  inner  region  where  no  solutions  are  available  from  the  difference 
scheme.  Therefore,  values  at  grid  points  in  the  inner  region  must  be 
obtained  by  some  other  method. 

Since  the  values  of  the  variables  are  available  at  the  inner/ 
outer  boundary;  we  use  a  method  based  on  extrapolations  from  the  boundary. 
The  position  of  the  boundary  at  each  x  directed  grid  line  is  determined 
by  a  linear  fit  to  the  two  nearest  neighbors  to  the  grid  line.  Values 
for  the  variable  for  the  two  components  of  velocity  are  obtained  in  a 
similar  manner.  The  depth  extrapolated  to  points  in  the  inner  region  is 
determined  by  (4.1.10).  Since  the  difference  scheme  (A1.6)  actually 
uses  mass  flux  values  instead  of  velocities,  the  velocities  at  inner 


region  grid  points  are  determined  by  matching  mass  flux  values  to  the 
boundary.  That  is,  the  grid  point  velocities  are  chosen  so  that  the 
mass  flux  at  the  grid  point  is  equal  to  that  at  the  interpolated  boundary 


APPENDIX  II  -  OTHER  COMPUTATIONAL  TECHNIQUES 


AII.l  -  Depth  Derivatives  at  the  Free  Boundary 

At  each  drift  point  on  the  free  boundary  the  depth  derivatives 
are  required  to  compute  the  acceleration  terms  in  (3.3.1).  These 
derivatives  are  calculated  from  the  surrounding  grid  point  depths. 

Since  the  drift  points  are  usually  not  on  the  grid  lines,  the  technique 
used  for  the  grid  point  derivatives  of  differencing  nearest  neighbors 
cannot  be  applied  directly.  The  method  used  must  be  based  either  on 
interpolation  or  extrapolation.  Two  approaches,  one  of  each  of  these 
types,  have  been  tested  for  use  in  this  model. 

The  first  method  tested  was  an  extrapolation  technique  based  on 
that  used  by  Kasahara  et  al.  (1965).  In  the  case  examined  by  Kasahara 
et  al.  the  cross  stream  flow  was  assumed  to  be  zero  and  there  was  no 
region  of  turbulent  entrainment.  Therefore,  the  frontal  layer  was  not 
divided  into  inner  and  outer  regions  as  is  done  here,  and  the  thickness 
of  the  frontal  layer  was  taken  to  be  zero  at  the  free  boundary.  Since 
the  free  boundary  was  then  entirely  outside  of  the  frontal  layer,  extrap¬ 
olation  was  required. 

The  extrapolation  procedure  used  for  this  model  involves  con¬ 
structing  a  line  which  is  locally  normal  to  the  front  at  the  boundary 
point  for  which  the  derivatives  are  being  computed.  The  method  of 
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constructing  the  normal  line  is  discussed  in  All. 2.  The  points  at  which 
the  normal  line  crosses  grid  lines  in  the  frontal  region  are  noted,  and 
from  these  a  base  point  is  selected.  The  base  point  is  the  point  closest 
to  the  front  which  has  at  least  three  of  its  four  nearest  neighbors  on 
the  same  grid  line  inside  the  frontal  region  (type  1  or  2  points  -  see 
All. 2).  Values  for  the  depth  and  the  x  and  y  derivatives  of  the  depth 
are  obtained  at  the  base  point  by  interpolation  from  the  three  neighbors. 
From  these  values  and  the  known  depth  at  the  boundary,  least  squares 
quadratic  fits  are  made.  The  x  and  y  components  of  the  depth  derivatives 
at  the  boundary  point  are  computed  from  the  fitted  curves. 

In  the  present  model  this  procedure  was  found  to  yield  poor 
quality  results.  The  calculated  derivatives  were  not  smoothly  varying 
along  the  front  and  were  occasionally  not  within  acceptable  error  bounds. 
The  reasons  for  the  success  of  this  algorithm  in  the  Kasahara  et  al. 
study  and  its  failure  here  are  largely  attributable  to  the  difference 
between  the  exponential  interface  used  here  and  the  i inear  interface 
assumed  for  the  earlier  study.  Large  variations  in  interfacial  slope 
with  distance  from  the  front  exist  for  an  exponential  profile,  not  all 
of  which  can  be  accommodated  by  a  polynomial  fit.  Since  the  various 
derivatives  are  extrapolated  from  points  at  differing  distances  from  the 
front,  divergent  slope  estimates  were  obtained.  The  extrapolation 
procedure  is  also  computationally  expensive.  To  surmount  these  problems, 
a  second  scheme  was  devised  which  uses  the  depths  which  are  available  in 
the  inner  region  in  an  interpolation. 
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A  boundary  point  may  be  viewed  as  being  inside  a  square  with  a 
grid  point  at  each  corner.  For  an  x  derivative,  two  pairs  of  grid 
points  are  formed  by  taking  points  on  the  same  x  directed  grid  line  to 
be  a  pair.  An  x  derivative  is  formed  for  each  pair  and  these  are  linearly 
interpolated  to  the  boundary  point.  The  same  procedure  is  used  for  y 
derivatives  with  the  pairs  being  formed  along  y  directed  grid  lines. 

This  procedure  gives  stable  and  consistent  values  and  is  used  for  all  of 
the  model  runs. 

All. 2  -  Grid  Point  Status  Determination 

As  discussed  above,  grid  points  on  one  side  of  the  inner/outer 
boundary  (in  the  inner  region.  Appendix  I)  are  treated  differently  than 
points  on  the  other  side  of  the  boundary  (in  the  outer  region.  Chapter 
2).  Since  the  inner/outer  boundary  can  move,  a  grid  point  initially  on 
one  side  of  the  boundary  can  be  on  the  other  side  at  some  later  time  by 
having  the  boundary  move  across  it.  The  model  must  therefore  keep  track 
of  the  position  of  the  boundary  relative  to  each  grid  point  (the  grid 
point  status)  at  each  time  step. 

|  The  technique  used  here  to  determine  the  status  of  a  grid  point 

j  follows  the  method  of  Neumann  and  Sproull  (1979).  This  approach  is 

i 

!  based  on  a  mapping  of  a  point  in  a  two  dimensional  Cartesian  space  to  a 

;  vector  in  a  three  dimensional  space.  The  point  (x,y)  may  be  represented 

j 

by  (wx,  wy,  w)  with  w  t  0.  Taking  two  points 

P  *  (P-|,  P2»  P3) 


v  =  (vr  v2,  v3) 


the  line  between  them,  from  p  to  v,  may  be  expressed  as  the  vector 


P2V3  ‘  P3V2 
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Y  * 

b 

' 

p3vl  “  P1V3 
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c 

_  Plv2  P2V1  _ 

(All. 2.1) 


The  position  of  a  third  point,  r  -  (r^ ,  r^),  with  respect  to 

y  may  be  determined  by  taking  the  dot  product  r  ♦  y.  If  r  •  y  =  0,  then 
r  is  on  the  line.  If  r  is  to  the  right  of  the  line  looking  from  p  to  v, 
then  r  •  v  will  be  negative;  if  r  is  to  the  left  then  r  •  y  will  be 
positive.  The  perpendicular  distance  from  r  to  y  is  given  by 


d  «  (r  •  y)/[r3(a2  +  b2)1/2] 


(All. 2. 2) 


In  the  model,  the  two  boundary  points  closest  to  the  grid  point 
of  interest  are  chosen  and  the  boundary  is  represented  by  the  straight 
line  between  them.  Since  the  boundary  point  spacing  is  very  small 
compared  to  the  scale  of  the  disturbance,  the  use  of  a  linear  fit  does 
not  introduce  a  significant  error. 

This  procedure  for  determining  the  distance  from  the  front  to  a 
grid  point  is  also  used  in  the  determination  of  the  depth  field  in  the 
inner  reqion  (see  Appendix  I).  The  errors  introduced  for  this  comDu- 
tation  by  the  straight  line  approximation  have  been  tested  and  shown  to 
be  small. 


To  keep  track  of  the  numerous  grid  points,  a  set  of  numeric 
designations  has  been  adopted.  The  five  types  of  grid  points  are: 


Type  1  -  a  point  in  the  outer  region  whose  four  nearest  neighbors 
are  also  in  the  outer  region. 

Type  2  -  a  point  in  the  outer  region  which  is  at  least  1/4  of  a 
grid  unit  away  from  the  boundary,  but  which  does  not  have  all  of  its 
four  nearest  neighbors  in  the  outer  region. 

Type  3  -  same  as  type  2  except  that  it  is  within  1/4  of  a  grid 
unit  of  the  boundary. 

Type  4  -  a  point  outside  the  outer  region,  at  least  one  of  whose 
neighbors  is  in  the  outer  region. 

Type  5  -  a  point  outside  the  outer  region,  all  of  whose  neighbors 
are  outside  the  outer  region. 

The  time  steps  used  in  the  model  are  small  enough  so  that  the 
position  of  the  front  changes  by  much  less  than  1/4  of  a  grid  unit  in  a 
cycle.  Therefore,  the  status  of  any  point  will  not  change  more  than  one 
classification  in  a  cycle. 

The  technique  of  mapping  a  point  from  two  to  three-space  is  also 
used  in  the  determination  of  the  derivatives  of  the  depth  field  at  the 
frontal  boundary  for  the  extrapolation  method  discussed  above.  Two  of 
the  manipulations  required  for  that  operation  are  locating  the  point  of 
intersection  of  two  lines,  and  constructing  a  line  normal  to  a  given 
line.  The  point  of  intersection  of  two  lines  y  and  u  is  found  by  taking 
the  cross  product  of  the  transpose  of  the  two  vectors 
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s=  Jx  yT  (All. 2. 3) 

Note  that,  when  remapped  into  two-space, 

-s  =  YT  x  J  (All. 2. 4) 


is  the  same  point  as  s. 

For  y  and  a  to  be  perpendicular,  we  must  have 

(YTJ)  •  a  =  0  (All. 2. 5) 

where 

1  !  0  0  j 

J  *  I  0  1  0  | 

|  0  0  0  j 

The  line  a  which  is  normal  to  y  and  passes  through  the  point 
**(*!•  V  v3>  is 

j  bv3 

r.  =  -av3 

(a v2  -  bv1 ) 

All. 3  -  Redistribution  of  Points  Along  the  Boundary 

The  quasi -Lagrangi an  drift  points  which  mark  the  position  of  the 
free  boundary  are  initially  distributed  approximately  evenly  along  the 
length  of  the  boundary.  Later,  when  the  front  has  evolved,  the  drift 


(All. 2. 6) 

I 

I 

I 


points  may  become  clustered  in  some  regions  and  widely  spaced  in  others. 

If  the  spacing  gets  to  be  too  uneven,  it  becomes  advisable  to  redistrib¬ 
ute  the  points  along  the  boundary.  This  type  of  problem  is  common  to 
lagrangian  models.  Fortunately,  in  this  case  the  rezoning  can  be  performed 
as  a  one  dimensional  problem,  whereas  most  Lagrangian  flow  models  must 
be  rezoned  in  two  dimensions. 

The  rezoning  is  applied  along  the  length  of  the  boundary  and 
results  in  a  new  set  of  y  axis  coordinates  for  the  drift  points.  These 
new  positions  are  then  used  with  a  series  of  fitted  curves  to  determine 
the  corresponding  x  axis  coordinates  for  the  drift  points.  Since  the 
drift  points  will  begin  to  shift  as  soon  as  the  model  execution  is 
resumed,  the  rezoning  of  the  y  axis  coordinates  does  not  need  to  be 
precise.  However,  once  the  new  y  coordinates  are  chosen,  the  x  coordi¬ 
nates  must  be  determined  accurately. 

A  segmented  front  is  established  as  a  series  of  straight  line 
segments  between  the  points  on  the  evolved  boundary.  The  average 
interval  spacing  is  computed  by  summing  the  lengths  of  all  of  the  line 
segments  and  dividing  by  the  number  of  intervals.  A  new  set  of  points 
is  then  established  along  the  segmented  front,  distributed  evenly  at  the 
average  interval  spacing.  The  y  axis  coordinates  of  these  points  are 
used  as  the  rezoned  y  coordinates. 

A  cubic  spline  fit  is  then  made  to  the  unmodified  evolved  boundary 
points.  Rezoned  x  coordinates  are  interpolated  at  the  rezoned  y  positions 
using  this  fit.  Thus,  the  rezoned  points  are  distributed  evenly  along 
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the  length  of  the  boundary  (by  the  choice  of  y  coordinates)  and  determine 
a  boundary  with  the  same  shape  as  the  original  (by  the  determination  of 
the  x  coordinates).  The  new  velocity  components  are  determined  by  the 
same  method  of  interpolation  from  a  cubic  spline  fit  as  the  x  coordinates. 
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area  into  the  ambient  region  and  the  inner  and  outer  frontal  regions.  The 
unperturbed  surface  front  lies  parallel  to  the  y  axis  (into  the  page). 
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Plan  view  of  the  initial  perturbed  position  of  the  inner/ 
outer  boundary  and  the  surface  front,  and  the  position 
of  the  inner/outer  boundary  at  a  later  time  after  frontal 
evolution. 
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Figure  3.  Model  grid  showing  the  positions  of  the  boundaries.  The 
dotted  line  through  the  solid  circles  (quasi-Lagrangian 
drift  points)  shows  the  position  of  the  free  boundary. 
The  dashed  lines  show  equivalent  pairs  of  grid  lines 
which  make  the  top  and  bottom  boundaries  periodic. 
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Figure  4 


A  section  of  the  free  boundary  showing 
components.  The  component  of  velocity 
to  the  front,  qj_,  is  controlled  by  the 
velocity,  qj_,  and  the  velocity  due  to 

qif  ^ 


the  velocity 
perpendicular 
particle 
the  mass  flux. 
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A  section  of  the  free  boundary  showing  the  mass  flux 
velocity  components.  The  grid  oriented  velocities, 
u,  and  9^,  are  components  of  the  flux  velocity,  qf, 
wnich  isequal  to  the  sum  of  v  and  v  ,  the  boundary 
perpendicular  and  parallel  velocities. 


Normalized  energy  vs.  time  (in  inertial 
coefficient:  a)  0.001,  b)  0.0015,  c)  0 


Depth  contours  for  the  standard  case  and  two  comparison 
cases.  Dashed  lines  show  free  boundary  position,  <$  = 

=  0.33.  Solid  lines  are  contours  at  5  =  0.4,  0.5,  0.6, 
0.7,  0.8,  and  0.9.  The  four  plots  are  for:  a)  initial 
state,  b)  grid  sensitivity  case,  c)  standard  case,  d)  time 
step  sensitivity  case,  all  cases  after  ten  inertial 


x  position  (grid  coordinates) 
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Time  (inertial  periods) 


Figure  9.  Frontal  position  vs.  time  for  unperturbed  front  with 
varying  cross  stream  Froude  numbers.  Fp  =  -0.02, 

F  =  0.0,  F  8  0.01  initially  at  x  =  3.5,  F  *  0.01 
initially  at  x  8  5.5.  Solid  symbols  and  connecting 
lines  show  the  corresponding  linear  fits  to  the 
frontal  positions. 


u  (dimensionless) 
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Figure  10.  Cross  stream  velocity  across  the  frontal  region  for  the 

unperturbed  front  after  ten  inertial  periods.  The 

solid  line  is  for  F  =  0.0,  the  dotted  line  is  for 

F  =  0.01  and  the  dashed  line  is  for  F  =  -0.02. 
r  r 


sequence  of  frontal  position,  large  amplitude  disturbance.  Time  advances 
e  right.  Curves  are  for  0,  2,  6,  and  10  inertial  periods. 


Spectral  Power  (arbitrary  units) 


Time  series  of  frontal  position  for  the  Sargasso  Sea 
short  wavelength  tests,  a)  F  =  -0.02,  b)  F  =  0.00, 
c)  F  =  0.01.  Time  advances  to  the  right  at  one 
inertial  period  intervals.  The  x  scale  for  each 
curve  has  been  expanded  by  a  factor  of  five. 


18.  Time  sequence  of  frontal  position  for  the  Sargasso  Sea 
long  wavelength  3  km  amplitude  test.  Time  advances  to 
the  right  at  one  inertial  period  intervals.  The  x  seal 
for  each  curve  has  been  expanded  by  a  factor  of  five. 


